{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 1D optimal transport\n",
    "\n",
    "\n",
    "This example illustrates the computation of EMD and Sinkhorn transport plans\n",
    "and their visualization.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Author: Remi Flamary <remi.flamary@unice.fr>\n",
    "#\n",
    "# License: MIT License\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pylab as pl\n",
    "import ot\n",
    "import ot.plot\n",
    "from ot.datasets import make_1D_gauss as gauss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate data\n",
    "-------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%% parameters\n",
    "\n",
    "n = 100  # nb bins\n",
    "\n",
    "# bin positions\n",
    "x = np.arange(n, dtype=np.float64)\n",
    "\n",
    "# Gaussian distributions\n",
    "a = gauss(n, m=20, s=5)  # m= mean, s= std\n",
    "b = gauss(n, m=60, s=10)\n",
    "\n",
    "# loss matrix\n",
    "M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))\n",
    "M /= M.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot distributions and loss matrix\n",
    "----------------------------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAADFCAYAAABzYARGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4U2X2wPHvoWW1yFZEsSCgiGylQEUERHBQdFBRB0cEFVdcBmVUEHQQEZfB0UH9CeIwMo47KCqDK4iAICrSYgsUZEepiiJiEdm6vL8/ToKltCVt09ws5/M8edIkNzcnN2nOfXdxzmGMMcZUVBWvAzDGGBMdLKEYY4wJCksoxhhjgsISijHGmKCwhGKMMSYoLKEYY4wJCksoxhhjgsISijHGmKCwhGKMMSYo4r0OoKjExETXrFkzr8MwxhgDpKen/+ScaxjItmGXUJo1a0ZaWprXYRhjjAFE5OtAt7UqL2OMMUFhCcUYY0xQBJRQRORcEVkrIhtEZHQxj1cXkRm+x5eKSDPf/VVF5HkRWSkia0Tk7uCGb4wxJlwcsQ1FROKAycDZQDawTERmO+dWF9rsOmCnc+4kERkIPAJcBlwKVHfOtReRWsBqEXnVObcl2G/EGFN+ubm5ZGdns2/fPq9DMR6pUaMGSUlJVK1atdz7CKRRvguwwTm3CUBEpgP9gcIJpT8wzvf3TGCSiAjggKNEJB6oCRwAdpU7WhOQOXNg9WoYPhyqWKWmCUB2dja1a9emWbNm6L+uiSXOOXbs2EF2djbNmzcv934C+bk5Htha6Ha2775it3HO5QE5QAM0ufwGfA98AzzmnPu56AuIyFARSRORtO3bt5f5TZjfPfMM/PGPcMcdMGgQ7N/vdUQmEuzbt48GDRpYMolRIkKDBg0qXEKt7PPXLkA+0BhoDtwpIi2KbuScm+qcS3XOpTZsGFB3Z1OEczBmDNx8M5x3Hjz4IMyYoX/n5HgdnYkElkxiWzA+/0CqvL4FmhS6neS7r7htsn3VW3WAHcAg4APnXC7wo4gsAVKBTRUN3Bzqtttg0iS49lr4178gPh6aNtXbPXvCkiWQkOB1lMaYaBZICWUZ0FJEmotINWAgMLvINrOBIb6/BwDznS5W/w1wFoCIHAV0Bb4KRuDmd19/DZMnw9Ch8OyzmkwArrwSXn8dVqyAV17xNkZjjuShhx6ibdu2JCcnk5KSwtKlS70O6TAJvrOy7777jgEDBpS43S+//MLTTz9d6r66desGwMKFCzn//PPLFMesWbNYvfr3ZuyxY8cyb968Mu2jMhwxofjaRIYBc4A1wGvOuSwRGS8iF/o2mwY0EJENwB2Av2vxZCBBRLLQxPScc25FsN9ErHvuOb2++24oWmrt3x/atoVp00IflzGB+uyzz3jnnXdYvnw5K1asYN68eTRp0uTITzyCvLy8IER3uMaNGzNz5swSHy8tofhj+vTTT8v9+kUTyvjx4+nTp0+59xcsAU294px7D3ivyH1jC/29D+0iXPR5u4u73wRPfr4mlD59oLgp0ETguuu0kX7VKmjXLuQhmgjz179CRkZw95mSAk88UfLj33//PYmJiVSvXh2AxMTEg4999NFHjBgxgry8PE499VSmTJlC9erVD07TlJiYSFpaGiNGjGDhwoWMGzeOjRs3smnTJpo2bcpLL73EqFGj+OCDD6hSpQo33HADt956K+np6dxxxx3s3r2bxMRE/vvf/3LccccdEtfmzZsZNGgQu3fvpn///gfv37JlC+effz6rVq0iKyuLa665hgMHDlBQUMAbb7zBvffey8aNG0lJSeHss8+mX79+3HvvvdSrV4+vvvqKdevWkZCQwO7duwHYtWsX/fr1Y8OGDfTu3Zunn36aKlWqHLLNzJkzeeeddxg6dCizZ8/m448/5sEHH+SNN97ggQce4Pzzz2fAgAGlHq8hQ4bw9ttvk5uby+uvv84pp5wSrI8YsJHyEe+jj+CbbzRplOTKK6FqVSulmPB1zjnnsHXrVk4++WRuueUWPv74Y0B7n1199dXMmDGDlStXkpeXx5QpU464v9WrVzNv3jxeffVVpk6dypYtW8jIyGDFihUMHjyY3Nxcbr31VmbOnEl6ejrXXnstf/vb3w7bz/Dhw7n55ptZuXLlYcnG75lnnmH48OFkZGSQlpZGUlISEyZM4MQTTyQjI4NHH30UgOXLl/Pkk0+ybt26w/bxxRdf8NRTT7F69Wo2btzIm2++WeJ769atGxdeeCGPPvooGRkZnHjiiQcfO9LxSkxMZPny5dx888089thjRzyOZRV2k0Oaspk2DerXh4suKnmbxESt+nrxRZgwAXwngcYUq7SSRGVJSEggPT2dxYsXs2DBAi677DImTJhAx44dad68OSeffDIAQ4YMYfLkyfz1r38tdX8XXnghNWvWBGDevHncdNNNxPsaF+vXr8+qVatYtWoVZ599NgD5+fnFJowlS5bwxhtvAHDllVcyatSow7Y5/fTTeeihh8jOzuaSSy6hZcuWxcbUpUuXEsd4dOnShRYttAPs5ZdfzieffFJqG01J1q5dW+rxuuSSSwDo3LlzqUmrvKyEEsF27IBZs7QEcqQkcf31uv3sot0pjAkTcXFx9OrVi/vvv59JkyYd/CEvSXx8PAUFBQCHjZ846qijSn2uc462bduSkZFBRkYGK1euZO7cucVue6TutIMGDWL27NnUrFmTP/7xj8yfP7/Y7UqLqehr+G8Xvj8Ysxj4qxTj4uIqpX3JEkoEe+klOHCg9Oouvz59oEkTq/Yy4Wnt2rWsX7/+4O2MjAxOOOEEWrVqxZYtW9iwYQMAL774ImeeeSagS12kp6cDlJp8zj77bP71r38d/AH9+eefadWqFdu3b+ezzz4DdOqZrKysw57bvXt3pk+fDsDLL79c7P43bdpEixYtuO222+jfvz8rVqygdu3a/PrrrwG//y+++ILNmzdTUFDAjBkz6NGjBwCNGjVizZo1FBQU8NZbbx3cvqT9l3a8QsESSoRyTrsIn3oqtG9/5O3j4uCaa2DuXG1zMSac7N69myFDhtCmTRuSk5NZvXo148aNo0aNGjz33HNceumltG/fnipVqnDTTTcBcN999zF8+HBSU1OJi4srcd/XX389TZs2JTk5mQ4dOvDKK69QrVo1Zs6cyahRo+jQoQMpKSnF9rp68sknmTx5Mu3bt+fbb4sOv1OvvfYa7dq1IyUlhVWrVnHVVVfRoEEDunfvTrt27Rg5cuQR3/+pp57KsGHDaN26Nc2bN+fiiy8GYMKECZx//vl069btkCq5gQMH8uijj9KxY0c2btx48P7SjlcoiA4XCR+pqanOFtg6spUrITkZpkyBQL8vW7ZA8+bajlJMVbCJYWvWrKF169Zeh2E8Vtz3QETSnXOpgTzfSigRyl9N269f4M9p1kzHpCxYUCkhGWNinCWUCLVgAZx4oraLlEXv3vDJJ5CbWzlxGWNilyWUCFRQAIsWQa9eZX9ur17w229gtYrGmGCzhBKBMjNh587yJRR/hw+r9jLGBJsllAi0cKFelyehJCZqrzD/PowxJlgsoUSgBQvgpJMgKal8z+/VS6ezP3AgqGEZY2KcJZQIk5+v7Se9e5d/H717w549sGxZ8OIypiJ27NhBSkoKKSkpHHvssRx//PEHbx+opDOf5cuX88EHHwS0bY8ePcjwzZjZt2/fUgctTpw4sdRR7ddccw1r164lLy+PunXrVijmt9566+BcYeHA5vKKMBkZugJjeaq7/Hr21OsFC6B796CEZUyFNGjQ4OAP9rhx40hISGDEiBEBPz8/P7/UwY3FWb58OatWreLcc88t0/PmzJlT6uMTJ07k2muvpUaNGsXG+ZxvvYnyTH1SNGb/AMhwYQklwlSk/cSvQQMdFLlwoS4bbMwhvJi/vhQXXHAB3333Hfv27eP222/n+uuvJy8vj8TERK6++mrmz5/Pv/71L7Zv387IkSNJSEigW7dubN26lVmzZrF7926GDRvG6tWryc3NPbh2yPjx49m7dy8LFy5kzJgxh0zGuGfPHoYMGcKqVato06bNISWOpKQkVq1aRVxcHH/+85/57rvvyM/PZ9y4cWzdupUff/yRM844g0aNGvHBBx8cFufIkSOZNGkS7XxrSdx222189NFHNG7cmOnTp9OgQQN69OjBpEmTSElJYdu2bfTo0YOVK1ceFvMvv/zCqlWreOKJJ9i8eTPXXnstO3bsoFGjRjz33HMkJSVxxRVX0KBBA5YtW8a2bdv45z//WWmJyKq8IsyCBXDyydC4ccX207u3tqPs3x+cuIypLM8//zzp6eksW7aMiRMnsnPnTgBycnLo2bMnK1asoEOHDtxyyy3MnTuXtLQ0tm3bdvD548eP59xzz+WLL75g/vz53HnnnYgIY8eOZfDgwWRkZBw2s++kSZOoV68ea9asYcyYMXz55ZeHxfXee+/RrFkzMjMzD85cfPvtt3PMMcewePHigysoFo7z9NNPP2QfOTk5dO/enaysLE4//XQeeOCBEo9DzZo1S435lltu4frrr2fFihVceumlh8zI/OOPP7JkyRJmzZrF3XffHeCRLzsroUSQvDxYvBgGDqz4vnr1giefhC++gDPOqPj+TBTxYv76Ujz++OPM9k2TnZ2dfXDhqmrVqh080169ejWtWrXihBNOAHQK+BdeeAGAuXPn8v777zNhwgRAZ+395ggT2i1atIi77roLgI4dO9K2bdvDtklOTmb06NGMHj2aCy64gO4l1B8XjrOo+Ph4Lr1U1yC84oorGDRoUKlxlWbp0qW88847AFx11VXce++9Bx+76KKLEBGSk5NLnJMsGKyEEkG+/BJ27apYg7xfz566mqONRzHhbN68eSxatIjPP/+czMxMkpOTD1Y/1axZ84hTy4NOVT9r1qyDU9V/8803B9cLqYjWrVuTlpZG27ZtGT16NA8//HCx2wUaJ/w+XX1pU/OXR/VC61tU5vyNllAiyKJFeh2M2ajr19d2FP8+jQlHOTk51K9fn5o1a5KVlcWyEromtmnThrVr17J161acc8yYMePgY3379uWpp546eNtffVXaFPM9e/bklVdeASAzM7PYqe2//fZbEhISuPLKK7nzzjtZvnz5EfdbVF5e3sGFrl555ZWD09YXnpq/8Nr1pe27a9euvPbaawC89NJL9PT3vgkhSygRJC0NmjaFElYiLbOuXXWfYTbhtDEH9evXjz179tCmTRvGjBnDaaedVux2tWrVYtKkSfTp04fU1FTq1q1LnTp1AJ3m/rfffqN9+/a0bduWcePGAXDWWWeRmZlJx44dD/nRBhg2bBg7duygdevWPPDAA3Ts2PGw18zMzOTUU08lJSWFhx9+mHvuuQeAoUOH0qdPH/r06XPE91enTh0WL15M27Zt+eSTTxjj6yUzcuRInnzySTp16nSwzehIMU+ePJmpU6eSnJzMjBkzePzxx4/4+sFm09dHkJYttVRxhIXsAvbss3DDDbB+vQ6UNLErGqav3717NwkJCTjnuPHGG2nfvj233nqr12FFFJu+Pkbs3AkbNkBqQB9rYPz7svxtosGUKVNISUmhTZs27N27lxtuuMHrkGKO9fKKEL7q2aAmlLZtdS36tLTg9BwzxksjR44MaHVEU3mshBIh/KWIzp2Dt8+qVXW8mZVQDFRu7x8T/oLx+VtCiRBpadCihfbOCqbUVEhP1zVWTOyqUaMGO3bssKQSo5xz7Nixo9jpYsrCqrwiRFoadOkS/P2mpsLkybBuHZxySvD3byJDUlIS2dnZbN++3etQjEdq1KhBUnmnMPexhBIBfvoJtmyBm28O/r79VWhpaZZQYlnVqlVp3ry512GYCGdVXhHAN74pqA3yfq1bQ82a1o5ijKm4gBKKiJwrImtFZIOIjC7m8eoiMsP3+FIRaVbosWQR+UxEskRkpYhUrJIuBvkTSqdOwd93fDx07Pj7axhjTHkdMaGISBwwGTgPaANcLiJtimx2HbDTOXcS8DjwiO+58cBLwE3OubZALyA3aNHHiLQ0HdRYxrV4Apaaqt2S8/MrZ//GmNgQSAmlC7DBObfJOXcAmA70L7JNf+B5398zgT+IznJ2DrDCOZcJ4Jzb4Zyzn60ySkurnOouv9RUXcHxq68q7zWMMdEvkIRyPLC10O1s333FbuOcywNygAbAyYATkTkislxE7iruBURkqIikiUia9TI51A8/wNatlZ9QwNpRjDEVU9mN8vFAD2Cw7/piEflD0Y2cc1Odc6nOudSGDRtWckiRpTIb5P1OPhkSEiyhGGMqJpCE8i3QpNDtJN99xW7jazepA+xASzOLnHM/Oef2AO8BldC0HL3S0nTdkmImOw2auDht8LeEYoypiEASyjKgpYg0F5FqwEBgdpFtZgNDfH8PAOY7HXI7B2gvIrV8ieZMYHVwQo8NaWnQqhXUrl25r5OaqsuI5+VV7usYY6LXEROKr01kGJoc1gCvOeeyRGS8iFzo22wa0EBENgB3AKN9z90JTESTUgaw3Dn3bvDfRvTKyKic7sJFdewI+/bB2rWV/1rGmOgU0Eh559x7aHVV4fvGFvp7H3BpCc99Ce06bMro55+1Qb5Dh8p/Lf9rZGbqLMTGGFNWNlI+jK1YodehSCinnALVqmlCMcaY8rCEEsYyMvQ6JaXyX6tqVS2Z+F/TGGPKyhJKGMvMhEaN9BIKHTpYCcUYU36WUMJYZmZoqrv8OnTQgZTbtoXuNY0x0cMSSpjKzYWsrNBUd/n5X8tKKcaY8rCEEqa++goOHAh9CQUsoRhjyscSSpjy/6iHMqHUqwdNmlhCMcaUjyWUMJWRAdWr6yj5UEpJsZ5expjysYQSpjIzoV07XQArlDp00NHy+/aF9nWNMZHPEkoYci70Pbz8OnTQhbayskL/2saYyGYJJQxt2wbbt4e2h5ef9fQyxpSXJZQw5G/D8KKE0qKFro1i7SjGmLKyhBKG/KWD5OTQv3aVKtC+vZVQjDFlZwklDGVmQrNmULeuN6+fkqIxOOfN6xtjIpMllDCUkeFNdZdfhw6QkwNff+1dDMaYyGMJJczs3Qvr1nmfUMCqvYwxZWMJJcysWgUFBd708PJr317XsbeEYowpC0soYcbLHl5+Rx0FLVtaTy9jTNlYQgkzmZlQu7Y2ynvJ1kYxxpSVJZQwk5mp3YWrePzJdOgAmzbBrl3exmGMiRyWUMJIQYEmFC/bT/z8MfjXtTfGmCOxhBJGtmyBX3/1tv3Ez3p6GWPKyhJKGPFiDZSSHH881K9vCcUYEzhLKGEkM1PbTtq18zoS7TZsDfPGmLKwhBJGMjK0u26tWl5Hojp0gJUrdTp7Y4w5EksoYSRcGuT9UlJ05P769V5HYoyJBJZQwsQvv2ijfDi0n/hZw7wxpiwsoYQJf/fccEoorVvrEsQ2Yt4YE4iAEoqInCsia0Vkg4iMLubx6iIyw/f4UhFpVuTxpiKyW0RGBCfs6BNOPbz8qlfXpGIlFGNMII6YUEQkDpgMnAe0AS4XkTZFNrsO2OmcOwl4HHikyOMTgfcrHm70ysyExERo3NjrSA7lXxvFGGOOJD6AbboAG5xzmwBEZDrQH1hdaJv+wDjf3zOBSSIizjknIhcBm4HfghZ1FPKvgSLidSSH6tABXnxR17hv2NDraExA9u2DtWth61a95OTomUqTJnDCCdC8efh90UxUCCShHA9sLXQ7GzitpG2cc3kikgM0EJF9wCjgbKDE6i4RGQoMBWjatGnAwUeLvDydtv4vf/E6ksMVbpjv08fbWEwp9u6FDz6A11+Ht9+G3btL3rZFC/jzn/WSkmLJxQRNZTfKjwMed86V8u0G59xU51yqcy61YQyeBq9bB/v3h1f7iZ/19Apz+/fD449DUhJccgnMnQuDBsGMGfD555Cdrcll3Tr46COYMkUHOz36KHTqBF27wqJFXr8LEyUCKaF8CzQpdDvJd19x22SLSDxQB9iBlmQGiMg/gLpAgYjsc85NqnDkUeTLL/U6nMag+DVsqLUl/hhNmHAOpk+He+7R/ubnnAN33glnnaVd84pq2VIvZ50FN90EP/0Er70GDz8MZ54JF1wAjzyivTCMKadASijLgJYi0lxEqgEDgdlFtpkNDPH9PQCY79QZzrlmzrlmwBPAw5ZMDpeeDjVqQJuiXR3CROfOGqMJEz//DBddpCWROnVgzhy9nHNO8cmkOImJcMstWnL5+9/h44/1jObppzVZGVMOR0wozrk8YBgwB1gDvOacyxKR8SJyoW+zaWibyQbgDuCwrsWmZOnpWrUU6G9BqHXurG28v/7qdSSGzz6Djh3h/fdh4kRYvlwTSXnVqgWjR+t0CH36aEPepZfqSFtjyiignzDn3HvAe0XuG1vo733ApUfYx7hyxBf1Cgq0OunKK72OpGSdO+tJa2Ym9OjhdTQxbMoUuO027a21ZAmcemrw9n3MMdqYP3Ei3H23Jqp337UqMFMmNlLeYxs26Jl/p05eR1Iyf2xW7eUR52DcOK2iOu88PQMJZjLxq1IFRoyAxYu119gZZ8AXXwT/dUzUsoTiMf+PdOfO3sZRmsaN4dhjLaF4Ij8fhg2D+++Ha6+FN9/UdpPK1LWrloDq1NFG/A8/rNzXM1HDEorH0tN1ipO2bb2OpHTWMO+B/Hy46iptKL/rLnj22dA1tLVooUnlpJOgXz+YNSs0r2simiUUjy1fDsnJULWq15GUrnNn+Oor+M3mOwgN5+Dmm+GVV7Rr7yOPhH4A4rHHwsKF+uFfdhnMmxfa1zcRxxKKh5zThBLO1V1+nTtrBwIb4BgCzmmJ5N//hr/9TRvJvVK3Lrz3HrRqpV2VP/vMu1hM2LOE4qGNG3WapUhJKGDVXiHx97/DY49pF94HHvA6GqhXT0fgH3cc/PGPdlZhSmQJxUOR0CDv17gxNGpkCaXSvfCClkquuAL+7//CZ56tY4/VKq+EBE0q3xadLMMYSyieSk+HatXCv0Ee9HetUydLKJXqk0/g+uu1Z9V//qPdeMPJCSdo9deuXXDhhdagZg4TZt/Y2LJ8ObRvr0klEnTuDKtXw549XkcShTZtgosv1qnlZ84M314a7dvrHGIZGdoDraDA64hMGLGE4pFIapD38zfM+5crNkGSk6OTM+bnwzvvaJtFOOvXD/75Tx0TM2aM19GYMGIJxSObN8POnZGXUMCqvYKqoEDn3Vm3Dt54Q2cEjgTDh8PQodqB4PXXvY7GhAlLKB6JpAZ5v6Qknc7eEkoQTZjw+xxavXt7HU3gROCpp+D00+Gaa2DNGq8jMmHAEopH0tO1mrxdO68jCZyIjZgPqrlztcpo0CCdXiXSVKumpZOjjtL2n127vI7IeMwSikf8s5BXr+51JGVz2mm6XLFNZV9BW7bA5ZdrF7+pU8One3BZHX+8rg65YYOWVGwtlZhmCcUDubk6iWu3bl5HUnbdumm1/9KlXkcSwQ4c0PXc8/K0Yfuoo7yOqGJ69dKpYd58U5cjNjHLEooHMjJg3z7o3t3rSMqua1c9mf70U68jiWCjRsGyZfDcc5HTCH8kd9yh1V6jRtnZRgyzhOIB/49xJJZQjj5ahyJYQimnWbPgiSd0oaxLLvE6muARgWnTtOfGZZdpF0YTcyyheODTT3XQcePGXkdSPt26aRuQjWkroy1btJ2hc2f4xz+8jib46tXT9pTvvrP2lBhlCSXEnNNlJiKxdOLXrZt26MnK8jqSCJKbCwMHahZ+7bXI640RqC5dNFn+7386F5mJKZZQQmzrVp1XL9ITCli1V5mMGaNtC9Om6eJV0Wz4cJ3r6667dLliEzMsoYSY/0c4Ehvk/Vq00JmHLaEEaO5cPWu/8UYYMMDraCqfiE5uecwx2p5ifcxjhiWUEPv0U+0l2r6915GUn4iWUiyhBGDbNp1apV272OpS26ABvPyyLvoTiYM2TblYQgmxJUt0cGColgavLN266Vi2H37wOpIwVlCgM/L++qvO0FuzptcRhVbPnjB2rK7x8uKLXkdjQsASSgjt3q2L3UVy+4mf/z3YirCleOwx+PBD7SYcCYveVIYxY+DMM+GWW2D9eq+jMZXMEkoILVumM5RHQ0Lp1EmncrJqrxIsXaorL156Kdxwg9fReCcuDl56Sb8sAwfC/v1eR2QqkSWUEPL/+J5+urdxBEONGpCaagmlWL/8oj+exx8f2fN0BUtSks4KsHw53H2319GYSmQJJYQ+/VRrPurW9TqS4OjWDdLS7KTzEM5pb66tW+HVV6Pnw66oCy+EW2/Vjgnvvut1NKaSWEIJkdxcWLwYevTwOpLg6dFDk8nnn3sdSRiZNk0HLj74YHQURYPpH/+AlBQYMkQHY5moE1BCEZFzRWStiGwQkdHFPF5dRGb4Hl8qIs18958tIukistJ3fVZww48cn3+unX3OOcfrSIKnd2/trTZnjteRhImVK/UsvE8fHdRnDlWjhvZ227dP14DJy/M6IhNkR0woIhIHTAbOA9oAl4tImyKbXQfsdM6dBDwOPOK7/yfgAudce2AIELN9B+fM0fbJP/zB60iC5+ij9STcEgrahe/Pf9YqrpdegipW+C9Wq1bwzDOwaBHcf7/X0ZggC+Rb3wXY4Jzb5Jw7AEwH+hfZpj/wvO/vmcAfREScc186577z3Z8F1BSRKJ3EqHRz5+rU73XqeB1JcPXtq22t27d7HYnH/vIXWLtWB/M1auR1NOHtiivg2mvhoYdg3jyvozFBFEhCOR7YWuh2tu++YrdxzuUBOUCDItv8CVjunDusCVdEhopImoikbY/CX6afftLG62iq7vLzv6cPP/Q2Dk/99786eG/sWDgrZmt1y+app6BNGxg8GL7/3utoTJCEpFwuIm3RarAbi3vcOTfVOZfqnEtt2LBhKEIKqXnztPNP375eRxJ8nTrpLBsxW+21YoUO2uvVC+691+toIketWtp5Yfdu7WJt7SlRIZCE8i3QpNDtJN99xW4jIvFAHWCH73YS8BZwlXNuY0UDjkRz5+pSEampXkcSfHFx2gY9d24MLn+RkwN/+pO2m7z6qh4ME7g2bXSczqJFcM89XkdjgiCQhLIMaCkizUWkGjAQmF1km9nZ9NAgAAAQpklEQVRoozvAAGC+c86JSF3gXWC0c25JsIKOJM7p2XufPtH7e9O3r86BuHKl15GEkHNw9dWwebOeaR97rNcRRabBg7WE9+ijuia9iWhHTCi+NpFhwBxgDfCacy5LRMaLyIW+zaYBDURkA3AH4O9aPAw4CRgrIhm+yzFBfxdhLCtLF7CLxuouP387SkxVez32mC7n++ij0TW4yAsTJ+rCXFdfDevWeR2NqQBxYVZPkZqa6tLS0rwOI2j++U8YMQK++QaaNDny9pGqXTs9SY+JTjvz5sG558LFF2vpJNanVgmGb77RBrlGjXTQVu3aXkdkfEQk3TkXUIW9dZavZHPnQuvW0Z1MQEtgixfDnj1eR1LJNm7U8SannKKLSFkyCY6mTTU5r12r68cUFHgdkSkHSyiVaO9ebW+M5uouv7594cABWLjQ60gq0a+/Qv/+mkRmz7az6GA76yyd6+t//4P77vM6GlMOllAq0bvv6iwT/fp5HUnlO+MM/X194w2vI6kkBQV65vzVV3omHe3rwntl2DAd9Pjgg/D6615HY8rIEkolevllbVfo3dvrSCpfzZrapDBzpibRqHPPPXrmPHFidM2fE25E4OmndSrrIUPgiy+8jsiUgSWUSrJzJ7z3no7ZitbuwkUNHgy7dun7jirPPAOPPAI33aSTP5rKVb06vPWWno1dcIF2zTYRwRJKJXnjDW1TGDTI60hC56yz4JhjtGQWNd55R+fp6tdPpwuxRvjQOOYYeP99XffhvPPg55+9jsgEwBJKJXnlFWjZMjpHx5ckPl5LZO++q4sWRrz0dLjsMl3DY/p0fYMmdFq10mrGzZu1M0RU1qVGF0soleDbb7W306BBsXdCO2iQLroV8YOeV6/WsSYNG2opJSHB64hi0xlnwPPPwyefaHft3FyvIzKlsIRSCaZP15k5Yqm6y69LFzjxRC2hRaxNm+Dss7VE8uGHcNxxXkcU2wYOhMmT4e234aqrID/f64hMCSyhVIKXX9aqrpNP9jqS0BPRRDp/vk45E3Gys7UX1759mkxatvQ6IgM639cjj+jZ2o032sDHMGUJJcjWrIEvv9QeT7Fq8GAtoc2Y4XUkZeRPJjt26MRk7dp5HZEp7K67YMwYmDZNx6tYUgk7llCC7JlntKbkssu8jsQ7rVrBqafqsYiY2olNm7S+/vvvtXdRLPWmiCTjx2timTIFrrsugr5gscESShD9+CP8+986oDrWq91HjtSJY996y+tIAvDVV9Czpw6imT8funf3OiJTEhGYMAHGjdOVMgcPtob6MGIJJYiefFKr3keN8joS711yibYhPfxwmC+8lZYGZ56pP0oLF1rJJBKI6Fxfjz6q9aoXX6wrPxrPWUIJkpwcmDQJBgzQKp9YFxcHo0dre1LYrpPy1ltaMqlVS2fxbN/e64hMWYwYofWq77+vn+O3RReSNaFmCSVInn5aa0zuvtvrSMLH4MGQlKSllLDinC6Q9ac/QXKyrr9hZwGR6cYbdZzQ+vVw2mmQkeF1RDHNEkoQ7Nmjs26fey507Oh1NOGjWjVtS1m8WC9h4bffdGXAkSO1OLlggS7qZCLXeefpwEcRXT0zogdBRTZLKEHw7LOwfbtOSGsOdf31kJgIDz3kdSTo6PcuXeDFF7VRd/p0nSbZRL4OHXRm4o4dtWh84402VYsHLKFU0Lffavtgr17a69QcqlYtLQzMmePhWinOwXPPaV/mn37SZTTvuw+q2Nc/qhx3nJY4R42CqVOha1cdGGZCxv6jKsA5uOEGnbtq6lSvowlft9+uy4XffLOW5EJq61adKfjaazWhfPkl9OkT4iBMyMTHa7fid97RgaopKfD3v0NenteRxQRLKBXwn/9oB5MJE2yGjtJUrarz++XkaFIJSTfiggLN8m3bwscfa5/u+fOhceMQvLjxXL9+kJWl66ncc4+WVqzBvtJZQimnb77RM+9evXQWCFO6du3g/vu12qvSp2T5+GMdT3LjjXq9ciXcdptVccWaRo10CdHXX9d/2E6dtErhhx+8jixq2X9YORw4oB2FnNNSiv1OBWbECO3Z+Ze/wJYtlfACWVnaFbhXL20reeUV+OgjW/891g0YAGvXwvDhOrr+pJO0l8iuXV5HFnXsp7CMDhzQZRkWLNAF/Jo39zqiyBEfr1VfBQXQu3cQk0p6ug7Nb9dOG9wffFB/QC6/PPYWpDHFq1dP+/ZnZenSomPGwAknaOcMWw0yaCyhlIE/mfzvf5pMrr7a64giT6tWOiv8L79UMKkcOKB1Z717a7XWggUwdqzu8G9/s+7Apngnn6z/wMuW6Xdn/Hho2hRuukk7bJgKsYQSoD17Dk0m1m5SfqmpMG+eJpVevWDDhgCf6JzOvTVyJDRpogsvbdmi62R8/bU20jRoUImRm6iRmqrLiq5cqf/YL7ygbSynnabTXlg7S7mIC7OZ+1JTU11aWprXYRxi3jwYOlSXtp40SdsATMWlp+vCiPv3wwMPaLv5Ycu2HzgAn36q3elef10/hPh47cVz8826A2vEMhW1c6cOeJ06VavFqlTRSUMvvhj69tVunDFafSoi6c65gGZNtYRSim3bdG6u//5XS8r//rfOQWeCJztbF+N7+209aZzy+D5SJV2TyOLFWpW1e7cmkT/8Qc8mL7oI6tf3OnQTjZzThPL661qlunat3t+smZ68dOumlxhKMEFPKCJyLvAkEAc865ybUOTx6sALQGdgB3CZc26L77G7geuAfOA251ypc896nVDy87Vd99//1h8553Q9n7FjoUYNz8KKPvn52pVz/Xrcqiy2zF7BriUraJ23kmro+hYFLU6iSt+z4ZxztCH16KM9DtrEnI0b9Qdhzhxd3iAnR+9v0EAHTSYn6yzVrVtrkqlfP+oSTVATiojEAeuAs4FsYBlwuXNudaFtbgGSnXM3ichA4GLn3GUi0gZ4FegCNAbmASc750pcZi2UCaWgQBfo27BBJ5xdskRPjHfsgIYNYcgQreqyQYtlsG+fVh/88ot23f3xR718/70WR7Zu1USyefOhCyMdeyy5bZLJlI68sP50ZnzTlZzqjTj1VF3vqnt3/Z9t2lQnnTQm5AoKdCqXzz6DpUshMxNWrYK9e3/fpm5d7frZpIlekpJ0PMwxx+ilfn3dpm5dXeMhAgQ7oZwOjHPO9fXdvhvAOff3QtvM8W3zmYjEA9uAhsDowtsW3q6k16tIQtm+9mcyH36X/HwOXnJz4UAuHNgP+/bDb7t1wtldu/T3LrfQjAzHHavJI7kDdO5UTH1+ZSrtcyj82JH+Lnxd9FJQ8Pt14UvhA5afr9NU5OXpwfNf9u/Xy759etm7V3sq7NkDv/76+2X//uLfQ5UqOtdSUpL+o514oh7sE0+ENm30n63QW1m6VGsdliyB5ct/zz1VqugujjtO/zfr14c6dbRTl/9Stape4uP1fzYuTp8ncvjFr6S/jSmNFOST8MNGjt62jto/bODoH9aTsH0ztXZmU2vHVqrv+aXE5+bWSCC3Rm1ya9Qmr0YCedVqkV+tJvlVa5JftQb5VatTULUG+fHVcHFVKfBdXJU4CuLicVXicBKn11XicFIFqlTBSRWcCPivq1al69Tryv8ey5BQAvnJPB7YWuh2NnBaSds45/JEJAdo4Lv/8yLPPb6YgIcCQwGaNm0aSNzF+in9a/q8cFW5n8823yVcplr3ioj+GletqsWBqlWhevVDL7VqQe3amghq1/79Ureu9vmvV0+nGfafmSUmBpyhRXSmjK5d9fbevZpU1q/Xgs3mzdoJZ/t2reLOydFt9u4N89UhTRSKA072XQ5Xi99oyHaO4Uca8QP12HnwcvS+XdTe9ytHs4sEdlOTvdTkV2rxA9XZT3X2U5P9VOMAVck9eB1PiRU8xdpDTahAQimLUJ6Dl8g5NxWYClpCKe9+TrqwDT99voH4eA5eqlWLoE5ApZ0aB3I67f+78LX/UvgUvfApu/9v/6m8/3YYqVnz92qv0jinncJyc38vYOXnH1oYK1xgK/y84v42puKO8l2alfmZBcBe3+UQ/lqG/HykIP/368L383uthAicUNG3EaBAEsq3QJNCt5N89xW3TbavyqsO2jgfyHODpmpCdRJPO7Gydm/CnMjvBShjopegJaPwa4MJ5FR0GdBSRJqLSDVgIDC7yDazgSG+vwcA8502zswGBopIdRFpDrQEvghO6MYYY8LJEUsovjaRYcAcNCX+xzmXJSLjgTTn3GxgGvCiiGwAfkaTDr7tXgNWA3nAX0rr4WWMMSZy2cBGY4wxJSpLL6/wan01xhgTsSyhGGOMCYqwq/ISke3A1xXcTSLwUxDCiWR2DJQdBzsGfnYcyncMTnDONQxkw7BLKMEgImmB1vlFKzsGyo6DHQM/Ow6VfwysyssYY0xQWEIxxhgTFNGaUKZ6HUAYsGOg7DjYMfCz41DJxyAq21CMMcaEXrSWUIwxxoSYJRRjjDFBEVUJRUTOFZG1IrJBREZ7HU+oiEgTEVkgIqtFJEtEhvvury8iH4rIet91Pa9jrWwiEiciX4rIO77bzUVkqe87McM3wWlUE5G6IjJTRL4SkTUicnqsfRdE5Hbf/8IqEXlVRGrEwndBRP4jIj+KyKpC9xX72Yv6P9/xWCEinSr6+lGTUHxLFU8GzgPaAJf7liCOBXnAnc65NkBX4C++9z4a+Mg51xL4yHc72g0H1hS6/QjwuHPuJGAnEJqVhrz1JPCBc+4UoAN6PGLmuyAixwO3AanOuXbopLYDiY3vwn+Bc4vcV9Jnfx46A3xLdIHDKRV98ahJKOi69Rucc5uccweA6UB/j2MKCefc98655b6/f0V/QI5H3//zvs2eBy7yJsLQEJEkoB/wrO+2AGcBM32bxMIxqAP0RGcAxzl3wDn3CzH2XUBnUq/pW5+pFvA9MfBdcM4tQmd8L6ykz74/8IJTnwN1ReS4irx+NCWU4pYqPmy54WgnIs2AjsBSoJFz7nvfQ9uARh6FFSpPAHehi92BLkP9i3Muz3c7Fr4TzYHtwHO+qr9nReQoYui74Jz7FngM+AZNJDlAOrH3XfAr6bMP+m9mNCWUmCciCcAbwF+dc7sKP+Zb8Cxq+4iLyPnAj865dK9j8Vg80AmY4pzrCPxGkeqtGPgu1EPPvpsDjdE1eItWA8Wkyv7soymhhHS54XAjIlXRZPKyc+5N390/+IuwvusfvYovBLoDF4rIFrS68yy0LaGur9oDYuM7kQ1kO+eW+m7PRBNMLH0X+gCbnXPbnXO5wJvo9yPWvgt+JX32Qf/NjKaEEshSxVHJ11YwDVjjnJtY6KHCSzMPAf4X6thCxTl3t3MuyTnXDP3s5zvnBgML0GWpIcqPAYBzbhuwVURa+e76A7piasx8F9Cqrq4iUsv3v+E/BjH1XSikpM9+NnCVr7dXVyCnUNVYuUTVSHkR+SNaj+5fqvghj0MKCRHpASwGVvJ7+8E9aDvKa0BTdEmAPzvnijbYRR0R6QWMcM6dLyIt0BJLfeBL4Arn3H4v46tsIpKCdkyoBmwCrkFPHmPmuyAi9wOXoT0gvwSuR9sHovq7ICKvAr3Qaep/AO4DZlHMZ+9LtpPQ6sA9wDXOuQotlxtVCcUYY4x3oqnKyxhjjIcsoRhjjAkKSyjGGGOCwhKKMcaYoLCEYowxJigsoRhjjAkKSyjGGGOC4v8B1ZcTbJpyYBQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 460.8x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcHHWd//HXp3tyEBIIEAyBAANyLSSAIcrNAkFBwI0/UeRQoiL3fSyni6DLobIQlgU0giy3CAgIcqywIAgkmhCW+wwEEo4EDHckycz398enerpn0pPpmenqb1f3+/l49KNTXTVVHzrMp9+p/lZ9LYSAiIjUXi52ASIizUoNWEQkEjVgEZFI1IBFRCJRAxYRiUQNWEQkEjVgkSZhZtub2Qux65AiNWBpema2n5lNN7OPzewtM7vbzLbr5z5fM7NdqlVjBccLZrbesrYJITwcQtiwj/t/zcwWmdmILq/PTI7d2pf9Njs1YGlqZnY8MBk4BxgJrAVcCkyMWVe1mVlLFXbzKrBvyT7HAkOqsN+mpQYsTcvMVgR+AhwRQvh9COGTEMLiEMIdIYR/TbYZZGaTzezN5DHZzAYl60aY2Z1m9r6Z/d3MHjaznJldgzfyO5JUfVKZY+9oZnPM7CQzm5ck76+b2e5m9mKyv9NKtv+SmT2WHOstM/svMxuYrHso2ez/kuN9u2T/J5vZ28CVhdeSn/l8coxxyfLqZjbfzHZcxlt2DXBAyfIk4Oo+vfkCqAFLc9saGAzcuoxtTge2AjYHNgO+BPwoWXcCMAdYFU/PpwEhhPBd4HXgayGEoSGEn3ez79WS468BnAH8GvgOsAWwPfBvZrZOsm0bcBwwIql7AnA4fsAdkm02S453Y8n+VwbWBg4uPXAI4RXgZOBaMxsCXAlcFUJ4cBnvxVRgBTP7JzPLA/sA1y5je+mBGrA0s1WAd0MIS5axzf7AT0II80II84GzgO8m6xYDo4C1k+T8cOjdzVUWA2eHEBYDv8Wb60UhhI9CCM8Az+JNnxDCjBDC1BDCkhDCa8CvgH/uYf/twI9DCJ+FEBZ2XRlC+DXwMjAt+e84vYKaCyn4y8BzwNwKfka6oQYszew9YEQP50dXB2aXLM9OXgP4Bd7A/sfMZpnZKb09fgihLflzoUG+U7J+ITAUwMw2SE53vG1mH+LnrDt9IVbG/BDCP3rY5tfAGODiEMJnFdR8DbAf8D10+qHf1IClmT0GfAZ8fRnbvIn/E75greQ1kqR6QghhXeBfgOPNbEKyXbVvM3gZ8DywfghhBfx0h/XwM8uswcyG4l9AXgGcaWYr91RECGE2/mXc7sDvK6hblkENWJpWCOED/NzrJckXYEPMbICZfdXMCudtbwB+ZGarJkOwziA572lme5rZemZmwAf4edr25OfeAdatYrnDgA+Bj81sI+CwLuv7cryLgOkhhB8CfwR+WeHPHQjsHEL4pJfHky7UgKWphRD+Azge/2JtPvAGcCRwW7LJvwPTgSeBp4DHk9cA1gfuAz7G0/SlIYQHknXn4o37fTM7sQqlnoj/0/8j/LTBjV3WnwlclRxv7552ZmYTgd0oNvLjgXFmtn9PPxtCeCWEML0XtUs3TDdkFxGJQwlYRCQSNWARkUjUgEVEIlEDFhGJpBo36JAGMGLEiNDa2hq7DJFMmTFjxrshhFX7+vNqwAJAa2sr06drZJFIb5jZ7J636p5OQYiIRKIGLNJMPvgAHn0UZs8GXQMQnRqwSKMLAW66CcaMgeHDYdttobUVPvc5OOUU+Oij2BU2LTVgkUb29tswYQLsvTfk83D22XD77XDppbDTTvCzn8GGG8K998autCnpSziRRjVnjjffOXPgkkvgkEO8CRccdhhMmwYHHQRf+xrccAPstVe8epuQErBII5o7F7bf3hPw//wPHH545+ZbsOWW8PDD8MUvekq+6aba19rE1IBFGs3ixd5M330X7r/fz/kuy4or+imIrbeGSZPg6adrU6eoAYs0nJNO8pEOV1wB48dX9jNDh8LNN3sz3msv+PDDdGsUQA1YpLHceSdMngxHH+0puDdWWw1uvBFeeQWOOCKd+qQTNWCRRvHxx36ud8wY+MUv+raPHXaA006Da6+FP/2puvXJUtSARRrFmWfCG2/Ar34FAwf2fT+nnQbrr+/N/B89zekp/aEGLNIInnzSTz0cdBBss03/9jV4sI8TfvllOPfc6tQnZakBizSCk06CFVaA886rzv522QW+/W0/lfHmm9XZpyxFDVgk6x54wIeRnXYarNzjzPKVO+ccWLIEfvKT6u1TOlEDFsmyEPx+DqNHw5FHVnff667rV89dfjm8+GJ19y2AGrBItt16K/z1r3DWWX7uttp+9CPf77/9W/X3LWrAIpkVAvz7v/uIhQMOSOcYI0fCscf6JcrPPZfOMZqYGrBIVt19N8ycCaeeCi0p3lfr2GNhueWq9wWfdFADFsmiQvpday34znfSPdaIEXDwwXDddTBrVrrHajJqwCJZ9Oc/w2OP+fCzAQPSP96JJ/rd1H7+8/SP1UTUgEWy6Oc/9xktfvCD2hxvjTX8Tmn//d8wb15tjtkE1IBFsubZZ/3875FH+rnZWjn+ePjsM79KTqpCDVgkayZP9qFhhx5a2+NutBHssYc34IULa3vsBqUGLJIl8+bB1Vf7sLNVV6398U84AebP9y/kpN/UgEWy5Je/9NMAxx0X5/g77gibbw4XXKBp7atADVgkKxYtgssug69+1U8HxGDmzf+553y6I+kXNWCRrLj5Zp9k86ij4tbx7W/76Y+LL45bRwNQAxbJiosv9suOd901bh2DBvmFGXfcAa++GreWjFMDFsmC6dNh6lQfeparg1/bww7zOjQkrV/q4G9SRHp08cU+c/H3vhe7ErfGGj578uWXw6efxq4ms9SARerdu+/6bMUHHOCzXtSLI4+E99+HG26IXUlmqQGL1LsrrvChZ4cfHruSzrbbDsaOhUsu0ZC0PlIDFqlnbW0+9GzHHWGTTWJX05mZfyjMnOnnp6XX1IBF6tldd8Hs2XDEEbErKe873/HTIv/1X7ErySQ1YJF6dsklsPrqMHFi7ErKGzrU75J2003wzjuxq8kcNWCRevXSSz7b8SGH1Oaev311+OGweLGPiJBeUQMWqVeXXeZTDR10UOxKlm2jjWDCBPjVr3wae6mYGrBIPfr0U7jySh9rO2pU7Gp6dsQR8MYbcOedsSvJFDVgkXp0/fU+xrZev3zr6mtfgzXX9HPWUjE1YJF6E4I3srFjfaxtFrS0+A3i77sPnn8+djWZoQYsUm8eeQSeeMKvNDOLXU3lfvhDGDhQQ9J6QQ1YpN5cfDEMHw777x+7kt753Odgn33gqqvgww9jV5MJasAi9WTuXLjlFjjwQFh++djV9N5RR8HHH/vsydIjNWCRevLLX0J7e/3d96FS48fDllv6aYj29tjV1D01YJF68Y9/+FjaPfeEddeNXU3fHXWUX0Ryzz2xK6l7asAi9eK663zG4WOPjV1J/3zrWz52+cILY1dS99SARepBCN6wNt0UdtopdjX9M3Cgj+C47z546qnY1dQ1NWCRenDfffDMMz7jcJaGnnXnkENgueXgootiV1LX1IBF6sGFF8LIkbDvvrErqY5VVvEZPK69FubNi11N3VIDFontqafg7rv9suNBg2JXUz3HHQeLFmn6+mVQAxaJ7ec/9zG/WbnvQ6U23NDvY3zJJT42WJaiBiwS0+zZPqnlQQfByivHrqb6Tj4ZFiyAX/86diV1SQ1YJKYLLvAv3Y4/PnYl6dhqK/jnf/b/zkWLYldTd9SARWKZN8+T4f77+60cG9XJJ8OcOXDNNbErqTtqwCKx/OIXPt38qafGriRdu+0GW2wBZ5/tUxdJBzVgkRjmzYNLL4X99vMvqxqZGZx5Jrz6qlJwF2rAIjGcf77f++FHP4pdSW3ssYdScBlqwCK19vbbPjRr330bP/0WFFLwrFm6VWUJNWCRWjvrLB8RcOaZsSuprT328FERZ57pk46KGrBITb34oo98OOQQWG+92NXUlplfdPLmm7pHREINWKSWTjvNb1JzxhmxK4lj++19BuXzzoP33otdTXRqwCK18tBDPt3QiSf6/GnN6txz/dLkH/84diXRqQGL1MKSJX6P3LXXhn/919jVxLXJJj7l0mWX+ezPTUwNWKQWLrnE73p24YUwZEjsauL76U/9lpVHHNHUc8epAYuk7c03/ZzvV74CX/967Grqw/Dh8LOfwaOP+jT2TUoNWCRNIfiIh0WLfKbgRpjtolomTfIv5Y47DubOjV1NFGrAImm69lq480445xxYf/3Y1dSXXA5+8xv/cDr4YP+wajJqwCJpmTMHjj4att3Wn2Vp663noyLuusubcZNRAxZJw+LFsM8+Pvrhyishn49dUf066iifCfqoo3xi0iaiBiyShtNPh0ce8avedOph2XI5uO46WGEF+Na3mmr6IjVgkWq76Sa/1++hh3oKlp6NGgXXXw/PPw8/+EHTDE1TAxappkcfhe9+F7bZxsf8SuV23tnvFXHTTf4viCbQErsAkYbxwgs+C/Caa8Ltt8PgwbEryp4TToBXXvF7Ray1Fhx2WOyKUqUGLFINzz/vCc7Mv9EfMSJ2RdlkBhdf7CNIDj/cv7w8+ODYVaVGpyBE+uvpp/1b/PZ2ePBBfenWXy0tcPPNfv/gQw7xC1galBqwSH/cfbef7zWDBx6AjTeOXVFjGDTI7xw3caIPTzvuOGhri11V1akBi/RFW5tf3bbnnn4xwV//Cv/0T7GraiyFJnzssTB5Muy+u0/n1EDUgEV6a9YsP+Vw+uk+bvWhh2D06NhVNaZ83keTTJni7/PYsXDbbbGrqho1YJFKffKJ39Vs4439PrZXXw033ABDh8aurPEddBDMmOEfdP/v//n54RdfjF1Vv6kBi/Tkgw98WNQ66/h9bL/xDXj2WR/vq7ub1c7GG8O0aXD++fDww758wAH+d5FRasAi5SxZAv/7v/C97/lVWqeeClts4ZcXX3+9TjnEMnCgjxV+8UU45hg/R7zJJrDDDj7d/YIFsSvsFQtNeAs4Wdr48ePD9OnTY5cRT1ubX0jxyCM+muHee+Hvf/fTC/vt58Ohxo2LXaV0NX++30XtiivgpZf8nPEOO8CECf48bhwsv3xqhzezGSGE8X3+eTVggSZowIsWeTqaPx/eeccH+s+eDS+/7I336afh009929VWg1128XONu+6a6i+wVEkI8Le/wa23wh//6NM/gZ8i2mADH6Gy/vrQ2upXKo4a5ROjrryy//328VSSGrBUxfhVVw3TJ05M/0Dd/f9WeL10fQhLP9rb/bmtrfhYssQfixfDZ5/5Y+FCf3z8MXz0kf+5nNGjYcMNYcwYT0tbbum/sDq3m23vvef35Zg5078wfeEF/7BdtGjpbVtaYNgwfwwZAsst50PgBg3yUx4DBvg2+XzxsfrqcMEFasBSHeMHDgzTazVVenfNrfB66Xqzzo9czp9Lfxnyef8lGTCg+Iuz3HL+GDbMTyMMH+6PVVeFkSO98Y4e7dtKc2hv93/9vP66P8+b5/8qWrDAP6Q/+sj/FbRwYfGDfNEi/2BfsqT4gd/e7rNb33tvvxuw7gUhbtNNoZFPQYjkcn7qYdSo2JV00CgIEZFI1IBFRCLROWABwMw+Al6IXUcPRgDvxi6iB1moEbJRZxZq3DCEMKyvP6xzwFLwQn++TKgFM5uuGqsjC3Vmpcb+/LxOQYiIRKIGLCISiRqwFEyJXUAFVGP1ZKHOhq9RX8KJiESiBCwiEokasIhIJGrATc7MdjOzF8zsZTM7JXY9BWa2ppk9YGbPmtkzZnZM8vrKZvYnM3speV6pDmrNm9lMM7szWV7HzKYl7+mNZjYwcn3DzexmM3vezJ4zs63r7X00s+OSv+enzewGMxtcD++jmf3GzOaZ2dMlr5V978z9Z1Lvk2bW4/1L1YCbmJnlgUuArwIbA/uaWb1M67sEOCGEsDGwFXBEUtspwP0hhPWB+5Pl2I4BnitZ/hlwYQhhPWABcGCUqoouAu4JIWwEbIbXWjfvo5mtARwNjA8hjAHywD7Ux/v438BuXV7r7r37KrB+8jgYuKzHvYcQ9GjSB7A1cG/J8qnAqbHr6qbW24Ev41frjUpeG4VfQBKzrtHJL+HOwJ2A4VdvtZR7jyPUtyLwKskX7iWv1837CKwBvAGsjF8cdiewa728j0Ar8HRP7x3wK2Dfctt191ACbm6F//EL5iSv1RUzawW+AEwDRoYQ3kpWvQ2MjFRWwWTgJKA9WV4FeD+EsCRZjv2ergPMB65MTpNcbmbLU0fvYwhhLnA+8DrwFvABMIP6eh9Ldffe9fr3SQ1Y6pqZDQVuAY4NIXxYui54zIg2jtLM9gTmhRBmxKqhAi3AOOCyEMIXgE/ocrqhDt7HlYCJ+IfF6sDyLP3P/rrU3/dODbi5zQXWLFkenbxWF8xsAN58rwsh/D55+R0zG5WsHwXMi1UfsC3wL2b2GvBb/DTERcBwMyvcZyX2ezoHmBNCmJYs34w35Hp6H3cBXg0hzA8hLAZ+j7+39fQ+luruvev175MacHP7G7B+8m3zQPyLjz9Ergnwb5SBK4DnQggXlKz6AzAp+fMk/NxwFCGEU0MIo0MIrfh7978hhP2BB4BvJpvFrvFt4A0z2zB5aQLwLHX0PuKnHrYysyHJ33uhxrp5H7vo7r37A3BAMhpiK+CDklMV5cU68a5HfTyA3YEXgVeA02PXU1LXdvg/7Z4Enkgeu+PnWO8HXgLuA1aOXWtS747Ancmf1wX+CrwM3AQMilzb5sD05L28DVip3t5H4CzgeeBp4BpgUD28j8AN+Hnpxfi/Jg7s7r3Dv4C9JPldegof1bHM/ffrUmQz2w3/J1ceuDyEcF6fdyYi0mT63ICTMaQv4kOD5uD/nN03hPBs9coTEWlc/bkh+5eAl0MIswDM7Lf4N5ndNuARI0aE1tbWfhxS+uv9933y1zXX7Pz6jBkz3g0hrFpY/nLuW73/ZO4627Hluix2s77L61bYTy7Xeb/JcnF9YRblLvvJ5TuWO7bN5zvvq+P1ws/6c8h1OXa+cw2h4/XulpPn5OdC3sqvT57bW7r+XOGZzsvJYdq7rO+6XNiuuL7zftpbOq9f6rlwnK7btYRu9hs6r0+eSV6nJWDJn63FR+rl88lzsjxggI80G5BvA2Bgiz8Pyheeff1yLYsBGJw8L5/3KeaXyyfLLZ8BMCTnrw/L/wOAocnzCrmFyfrPkmV/fVjHs+9nmIVk2d+EobnBAORWe6mb6bz7rj8NuNyYty27bmRmB+NXhbDWWmsxXTPvRnXSSXDxxUtPgGxms+NUJNK8Up+SKIQwheSemePHj9e9LyNbvBgGDEhp54XTWYV0GZJrE5KEGtqTJJTrsr69c4ItnBaz9mR9IWUmy4XUaYVLH3Jd9kNb8pzvSHTWlrxWSMIFbe2dFi0ZGBTo/HohCRdqKlzDZHRdLugoLllPl/XJ2uQyg/alfhMLW3b+yVyy3N7NcleFd6Q92S6XbNdedusydXVTT3G/5Y8byv7Zf6qNrqrchvq7uyQJ0164BiRJ0P3cbdlD9eNn63oMqZS3aBEMjHprGBEp6M9nRccYUrzx7gPsV5WqJDWffgrLL5/yQeoqCeeTbZMSlISTZSXhHnVJwmkk4D6XGEJYYmZHAvfi/5f/JoTwTNUqk1R8+CEMHRq7ChGBfn5GhBDuAu6qUi1SA++9B6usUqOD1UUSLp4P9m2TEpSEk+X0k3DXL34ym4RToEuRm8zcubDaarGrEBGowSgIqR9tbfDaa/DNb/a4aXXFTMJlRkb4tkkJSsLJspJwDErATeTpp2HJEhgzJnYlIgJ1+Zkgafnzn/15hx0iFRAjCS9jjLBvm5SgJJwsKwnXkhJwE7n2Wth446UvQxaROOros0DSNHUq/O1vfhlydDVMwpVcLefbJiUoCSfL1UvClYwR7rzcPElYCbgJLFoEhx0GI0fCAQfErkZECurgM0DSFAKcfjo88QTcdhusUOnlPGbFpJqWGiTh3tw3ApSE00jCvblarvNy4ydhJeAGFoLf/ez88z0BT5wYuyIRKaUE3KA++ACOPhquvhqOOAL+8z/7sJOOZJrhJNyHO6j59kkJSsLJcn+ScO/vG9F5uXGTsBJwA7rtNh/tcO21cMYZ/sVbTn/TInVHCbiBPPYYnHsu3HEHbLqpN+IvfrGPO7NcSRLNcBLux72EffukBCXhZLkvSbjvd1DrvNx4SVi5KOPa2uDWW2HbbWGbbeAvf4FzzvEZL/rcfEWkJpSAM+q55+DGG+G66+Dll6G11c/zfv/7VbzdZGGutQwn4WrMquHbJyUoCSfLlSfhatxLuPNy4yRhNeAMefFFb7q/+53f18HMLys++2z4xjegRX+bIpmiX9k69vHH8PDDcN998Kc/wVNPedPdbjv/Ym2vvWDUqHSObTnrSJpZTsLVnF/Ot09KUBJOlntOwtWcVaPzcvaTsBpwHVm8GKZN84Z7//1++fCSJT6H27bbwuTJfivJNdaIXamIVIMacETz5nnDLTweeww++cQD3xZbwAknwC67ePNdbrna11dIlplOwinMtOzbJyUoCSfLTZCEU1D/FTaIzz7zy4GnTvVmO3UqvPqqr8vnYexYmDQJJkyAHXeElVeOWq6I1IAacAo+/hiefNIb7syZ/vzkk35THPBTCFtt5ZcHb7mlp93UZyrurZJxwJlOwqnMtAxKwqVHVxLuq/qtLCPeeadzo505E156qdgrVl4ZvvAFOOYYb7ZbbgmjR8etWUTqgxpwhRYuhGef9ZEIpY+33y5u09rqzXb//WHzzf3Po0cXA1um5KyY+rKchFOZaRmUhLvsv6MaJeHeqL+KImtrg1mzlm60L7/c8bvJ4MF+r4Vddy022s02g+HD49YuItnS1A143rylG+0zz8Cnn/p6M/j85/0Lsn328eexY2G99ZYOOA2pcC41w0k4lZmWS/ejJNx5/x3VFJNwf+4l3HmfjZeE66eSFH36qTfWrs123rziNquu6s31oIOKjXaTTerwyzERaRgN14D//nf/Iuzxx4vPL75YDEuDB/u07HvsUWy0Y8f6dD1SZFa8iizLSTiVmZZBSbgXSbgWMy2XVrz0cv0m4fgV9FEI8OabSzfb118vbrPmmjBuXOfTB5//fJOcPhCRupeZBtze7qcRHnrI74/w0EPw1lvF9RtsAFtv7bM/fOEL/hgxIl69mZfLFVNXlpNwCjMt+/rkmErCnZa7ypHOTMvLOmaWknDdNuDFi2HGjGKzfeQRWLDA162xhl8tttVWnnA32wyGDYtarohIr9VVAw7Bbyh+9dV+y8UPP/TXN9jAb7e4/fZ++8XW1oyOrc0Ss470l+kknMJMy6Ak3Jsk3Jcxwp1q6qaWRkjCPR7RzNYErgZG4v8NU0IIF5nZysCNQCvwGrB3CGFBX4qYNQuuucYb76xZPvJgr73ga1/zWy+utlpf9ioiUt8qaflLgBNCCI+b2TBghpn9CfgecH8I4TwzOwU4BTi5twX8x3/AiSd6WJkwAc4809Ouhn9FVjKjcKaTcBozLYOSsJJwVfR4pBDCW8BbyZ8/MrPngDWAicCOyWZXAQ/SywY8ZYo33732ggsv9FELIiLNolet3sxagS8A04CRSXMGeBs/RVGxV1+FQw+F8ePh+uv9puNSP8ysI91lOgmnMdMyKAkrCVdFrudNnJkNBW4Bjg0hfFi6Lvj/nWV/i8zsYDObbmbT58+f3/H66NGw004+fveuu/pWvIhIllXU4s1sAN58rwsh/D55+R0zGxVCeMvMRgHzyv1sCGEKMAVg/PjxHU16wAC47Tb48pf9FMSuu/oNySdO9KvVJLKcdaQ5JWGUhPuRhNOYabnzdtlNwj0mYDMz4ArguRDCBSWr/gBMSv48Cbi9twcfNgzuuQdOPtnvzbDPPj7i4eCDfdxve09/EyIiGVZJa98W+C7wlJk9kbx2GnAe8DszOxCYDezdlwKGD4dzzoGf/hQefBCuugquuw5+/WtYaSUfhrbDDj4GeNw4T85SA5brSG9KwigJ9yMJp3Ev4U41dVNL9ZNw9VUyCuIvLP33UTChWoXk8z4MbcIEuPRSuP12eOABvxLujjt8myFD/HLjwgUZX/qShquJSHbV1ZVwBUOH+qwS++/vy2+/7VfIFe4DcdZZxSGeG27o930YN654DwhNaFkFpTNiKAkrCRfW9iEJpzmrRqeauqmlnpNwXTbgrlZbDb75TX8AvP8+PPoo/O1vPoriL3+BG24obr/22sWGPG6cz1qx+uq6fFlE6ksmGnBXw4fD7rv7o+Ddd70ZFx6PP+6jLApBZ6WVirek3HRTfx4zRjfx6VYuT8dnv5KwknDnvfYqCddifrlONXVTS3+TcBoy2YDLGTHCh7R9+cvF1z76CP7v//zx1FM+NfzVV/vrBa2tnW/MPnas3/yn6b/syxmFZqBGrEbcv0bc+4s1ikeup0ZcfQ3TgMsZNsxHUWy3XfG1EGD27KWnJ7r7bliS/MUNHAgbbbR0Y87sDMciUpcaugGXY+apt7XV77ZW8Nln8MILnpILTfnPf/YhcQXDh/tpi66nMlZYodb/FenzS5ELS0rCSsL9ScJ9v2y5eOT4STgNTdeAuzNokDfUTTft/PqCBfD0053T8nXXFe9VDD5LcmF6+sLzaqspLYvIsqkB92CllXzc8fbbF18LAd54w9PyE0/44/HH4eabi9t87nPFhrz55rDFFt6oM9OU8/mOtKUkjJJwv5Jw/2/gUzxyYyVhNeA+MIO11vLHnnsWX//gA2/KM2d6U545Ey64wKdXAh+fvOWW/thqK7+QZKWV4vw3iEh8asBVtOKKS6flRYvg2Wdh+nSYNg2mTvX7XxQC1gYbeDPeckvYZhs/BZLLld9/TZl1pCslYZSEO36geZNwGtSAUzZwYPE0xA9/6K99+KE35KlTvSnfc48PjwMfTrfTTrDLLn5Z9rrrZui0hYj0ihpwBCusADvv7A8oDo176CG4/35/3HSTr2ttLd4jY7fdanjKorTrKwkrCffPc1UpAAAO+klEQVQjCac55b2vz24SVgOuA6VD4w44wH/XX3ih2IxvuQWuuMIvDvnKV2Dvvf2+ySuuGLtyEekPNeA6ZOYXgmy0ERxxBLS1+X0vbrkFfvc7+OMf/dTGbrt5M95rrxRuYp/PLZWelIRREu5DEk5zeqPikbOZhOvh6x7pQT7vX9T94hfw2mvw2GPemGfMgO98x28+9NOfwnvvxa5URHpDDThjzLwZX3ABvP463HefjzE+4wyfVfrII/18cr/lcp6G8jk/aOkjn/dxwmaYmae6nPkNfHL54rLl/JEsF7fP+aOwv2S5Y31HDV32U3gPclZMmP5Cp/Ud+01TCJ2Tdmgvpl08CXck9dL17cEfHbsJnobb2/1R2G+yXFyfPLrup70tefhyx/Ztbf4o7K/waGv3R7J/aw9Ye8nxC+uT7a293dNwW4C2csvJo609eQQsWddpfXsgt8QfxZ8pPEgeyXK7J/9cWyBXsr7rcmG74np/FPaTW+Ipt7j/Lo/Ccbput8T80WW/aVADzrBczr+cu+uu4pROU6bAxhvD5Mn++yci9UsNuEGMGQO/+Q289BLsuCMcd5yPK37mmb7tL+RKkmqWk3DalIR7TMKlKTjLSTgNasANZu214c474frrYdYsb8KPPBK7KhEpR6MgGpAZ7Luv34Zzl1186Nof/uCnKyqWy3V8Q25dP6ezNDqiFiMjSvev0RGJ4uiIvt1Bzfdauhx/dET1KQE3sDXX9Is71l3Xp3OaOzd2RSJSSgm4wY0cCbfe6veYOOggH0Nc0anRfPGzOdNJuJZjhEv3ryScaKdv940o3bJeknD1KQE3gfXWg3PO8Vk/HnggdjUiUqAG3CQOPdRvh3nppRX+gFlx9EMyOiHkLHujI2KMEQaNjigdAVHBGOEsjI5Igxpwkxg8GCZNgttvh08+iV2NiIAacFP5yld84tGpU3veNpQm0Awn4ZIFJeFISbg3V8vVcxJOgxpwE9l6a3/+61/j1iEiTqMgmsiKK/qoiFdeqWDjvBW/0S58E5/v/HmdhdERdXEHtdL9N+HoiHRmWl76J2s1OqKalICbTGur38RHROJTAm4yI0bAW2/1vF3I5ZYe25nFJFxP9xIu3X8TJeF05pcr3TK7SVgJuMmstBIsWBC7ChGBXiRgM8sD04G5IYQ9zWwd4LfAKsAM4LshhEXplCnVMmQILFxYwYal54CznITrcVaN0v03QRJOd6bl0i2zl4R7k4CPAZ4rWf4ZcGEIYT1gAXBgNQuTdAweXGEDFpHUVZSAzWw0sAdwNnC8+SDMnYH9kk2uAs4ELkuhRqmiAQNg8eKetws5K8kVGU7C9Ty/XOn+GzgJpzHTMtQ+Caeh0gQ8GTiJ4ju0CvB+CKFwgd4cYI1yP2hmB5vZdDObPn/+/H4VK/3X0uIXY4hIfD0mYDPbE5gXQphhZjv29gAhhCnAFIDx48enHCekJ7lcMTAuS8jnoCPJJq9lMAlnYqbl0v03YhJOYaZlqH0STkMlpyC2Bf7FzHYHBgMrABcBw82sJUnBowHdbTYDKm3AIpK+HhtwCOFU4FSAJAGfGELY38xuAr6Jj4SYBNyeYp1SJWaVhbmQN0ozAGQzCVdlVg1QElYSTkV/9n0y/oXcy/g54SuqU5KkqdIGLCLp69WVcCGEB4EHkz/PAr5U/ZIkTRXfACxvJWe+MpyEqzm/HCgJ9yEJV3N+OV9Pl/XJ2pSTcBp0JZyISCS6F0STqTQBdx4HXJDBJJzGTMugJNyLJJzGTMu+ni7rk7UZSsJKwCIikSgBS1mhZDaJLCfhVGZaLlmvJJyUoSTcJ0rAIiKRKAFLWe0tttRMsFlMwtWYVQOUhPuVhKswq0bn5YLaJuE0KAGLiESiBCxlhZx1JIZMJ+Eqzi8HSsLNnITToAQsIhKJErCU5feCcJlOwinMtAxKwr1KwinMtNx5uSDdJJwGJWARkUiUgKWskIeu2SGbSTiFmZZBSbgXSTiNmZZ9+9om4TQoAYuIRKIELGX5OeDyCSBbSTiFmZZL96MknKzvPgmnMdOyr691Eq4+JWARkUiUgKUsTy3LPheWhSScykzLoCTcqyRc/ZmWS2uqXRKuPiVgEZFIlIClrPa8lYx/zG4STmOmZVAS7l0STmF+OYiQhKtPCVhEJBIlYCkr5MvdDSqDSTiFmZZLa1ISriAJpznTMmQ6CSsBi4hEogQsZZWeA1YSVhLuVxJOYaZlqH0SToMSsIhIJErAUla5c8CZTMJpzLQMSsK9SMJpzC8HtU/CaVACFhGJRAlYygolH81KwigJ9yMJpznTMmQ7CSsBi4hEogQsZYX80q9lMQmnMdOyHzM5hpKwPy8rCacxvxzUPAmnQQlYRCQSJWApqz3f/adzlpJwKjMtg5Jwb5JwmjMtQ82ScBqUgEVEIqkoAZvZcOByYAweUH4AvADcCLQCrwF7hxAWpFKl1FzIG+1JFlUSBiXhvifhVGZahghJuPoqTcAXAfeEEDYCNgOeA04B7g8hrA/cnyyLiEiFekzAZrYisAPwPYAQwiJgkZlNBHZMNrsKeBA4OY0ipfbaWyCXZIBMJ+FUZloGJeHKk3AaMy1DjCRcfZUk4HWA+cCVZjbTzC43s+WBkSGEt5Jt3gZGlvthMzvYzKab2fT58+dXp2oRkQZQyTngFmAccFQIYZqZXUSX0w0hhGBmZT+KQwhTgCkA48ePT/njWqrF7wXhspyE05lpGZSEe5GEU5hp2bdPSqhREk5DJXueA8wJIUxLlm/GG/I7ZjYKIHmel06JIiKNqccEHEJ428zeMLMNQwgvABOAZ5PHJOC85Pn2VCuVmiq9Ei7LSTiNmZZBSbhXSTiVmZah1kk4DZVeiHEUcJ2ZDQRmAd/H/+//nZkdCMwG9k6nRBGRxlRRAw4hPAGML7NqQnXLkXpR/l4QLktJOI2ZlkFJuFdJOJWZlqHWSTgNuhJORCQS3QtCygrL+GjOUhJOZX45UBIuqCAJpzLTcul+apSE06AELCISiRKwlNXe0vNssFlIwlWdVQOUhLuzjCScykzLUPMknAYlYBGRSJSApaxOV8JlOgn34r4RoCTcX+WScAozLfv65Ji1SsIpUAIWEYlECVjKCi0BOhKty2YSrvy+EaAkXDWlSTiFmZahMZKwErCISCRKwFKWXwnXOXlmMQmnMdOyb68kXJEQenffCKjbJJwGJWARkUiUgKWszrMiZzcJV2VWDVAS7o80ZlqGmifhNCgBi4hEogQsZYV8KEmgHa8mz0rCdGyvJFwxJeGlKAGLiESiBCxl+Thgl+UknMpMy6Ak3B9ZTcIpUAIWEYlECVjKKk3ABVlMwmnMtAxKwlWhJKwELCISixKwlJcPdJd9lISVhKsqK0k4BUrAIiKRKAFLeSXngLOdhPtzL+FidUrCSsJpUAOWsqzMKYhsNuLqTW/Ueb0acWrqtBGnQacgREQiUQKWsqylncLnc7aTcN9vZVn6U0rCSsJpUAIWEYlECVjKyufbSz73s5uE05zyvvN6JeHU1E0Srj4lYBGRSJSApax8SzGJZTkJpzvlfbE6JeEmSMIpUAIWEYmkogRsZscBP8QDwFPA94FRwG+BVYAZwHdDCItSqlNqbMCAJXT93yOLSTjdKe9Lj64k7MsNnIRT0GMCNrM1gKOB8SGEMfjf9D7Az4ALQwjrAQuAA9MrU0Sk8VR6DrgFWM7MFgNDgLeAnYH9kvVXAWcCl1W7QIljQL507GN2k3CaU96X/pSScOMn4TT0uOcQwlzgfOB1vPF+gJ9yeD+EUPjfcg6wRrmfN7ODzWy6mU2fP39+daoWEWkAPSZgM1sJmAisA7wP3ATsVukBQghTgCkA48ePT/HjUKppYEu5q3+UhJWE+5CE00zBXmTxWJ2OXd0knIZKsvUuwKshhPkhhMXA74FtgeFmVviNHA3MTalGEZGGVMk54NeBrcxsCLAQmABMBx4AvomPhJgE3J5WkVJ7g/LLuv5dSVhJuBdJuBbng0v3n1YSTkEl54CnATcDj+ND0HL4KYWTgePN7GV8KNoVqVUpItKAKhoFEUL4MfDjLi/PAr5U9YqkLgzK9xA7gSwk4d7cN6J0WUm4ikm4liMjSvdf7SScAl0JJyISie4FIWUt17K4F1vXbxJOc6JPUBKuKAnHGCNcuv8qJeE0KAGLiESiBCxlDe5VAi6ovyRciynvQUl4WUk46tVypfvvZxJOgxKwiEgkSsBS1vL5/tzYrn6ScF/uoObrlYT9uMlx+pGE6+K+EaX772sSToESsIhIJErAUtZy+b6cA+5KSVhJuDgOOPNJOAVKwCIikSgBS1nLt3xWxb0pCTd1Eq7HewmX7r/SJJwCJWARkUiUgKWsIblFKfzfoSTclEm4nmfVKN1/D0k4DUrAIiKRKAFLWcPy/yguZDgJpzHTsq9XEvbjJsdZRhLOxPxypfvvJgmnQQlYRCQSJWApa2hpAi7IYBLuz72EO++zu2MqCftxk+OUScKZmmm5dP9dk3AKlIBFRCJRApayVsgt7H5lppJw/2fV6LzP7o6pJOzHTY5TmoTTmGnZd0SquibhFCgBi4hEogQsZQ3JVXAlXCaScPXmlyvuc1nHVBL24ybHyeXSmWkZap+EU6AELCISiRKwlLVCrswoiO7UcRLuz9VynbZTEvbt+pCEU5lpuWR9zZJwCpSARUQiUQKWsob1JgEX1GESTmOmZV9WEobKknAqMy37C53WZzEJKwGLiESiBCxlDcv1Y0YMJeEyx23iJJzCTMu+mP0krAQsIhKJErCUNcwC9CcFg5Jw2eM2YRJOYaZlaIwkrAQsIhKJErCUNSzXAu1J3MpwEk5zVo1O2ykJ+3Zlk3AKMy1DQyRhJWARkUiUgKWsobnBQDIWWElYSbhfSTiFmZZL95PhJKwE3GS22gqOOSZ2FSICYKGGnwZmNh+YXbMDSm+sHUJYNXYRIs2kpg1YRESKdApCRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUjUgEVEIlEDFhGJRA1YRCQSNWARkUj+P4VunBEI1yzqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% plot the distributions\n",
    "\n",
    "pl.figure(1, figsize=(6.4, 3))\n",
    "pl.plot(x, a, 'b', label='Source distribution')\n",
    "pl.plot(x, b, 'r', label='Target distribution')\n",
    "pl.legend()\n",
    "\n",
    "#%% plot distributions and loss matrix\n",
    "\n",
    "pl.figure(2, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, M, 'Cost matrix M')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve EMD\n",
    "---------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VdW5x/Hvm4RBJgFBQECj4nAVUDEtiENVrFqHi7dYxeoVby04D6h1bK1tr9apiperVipaB7SKVmmpw1XUakXQUBwYZBBEQTGACDIIJFn3j3eniRDIdJJ1ht/nefZzcvbZ55w3B/LLytprr2UhBEREpOnlxS5ARCRXKYBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsEiWMLNDzWxO7Dqk9hTAIgkzO8vMPjCzdWa21MzuNbP2yWO/N7M1ybbRzDZVuf98E9QWzKzXto4JIbwRQtirAe8x1MymmtlaMytJvj7fzCx53MzsFjNbkWy3VDwm9aMAFgHM7HLgFuBnwPbAAGAX4CUzax5CODeE0CaE0Aa4CXii4n4I4QfxKndmVtDA518O3AXcBnQFugDnAgcDzZPDRgAnAfsBfYETgXMa8r65TgEsOc/M2gG/Ai4KIbwQQtgUQvgYOAUoBM6ox2sebmaLzezKpDX5uZmdZGbHmdlcM/vSzK6tcvx3zewtM/sqOfZ/zax58tjryWHvJS3uU6u8/lVmthR4sGJf8pzdk/fol9zfycyWmdnh1dS6PfBr4PwQwlMhhK+Dmx5COD2EsCE5dBjwuxDC4hDCEuB3wFl1/WykkgJYBAYCLYE/V90ZQlgDPAd8v56v2zV53e7A9cAf8DA/EDgU+IWZ7ZocWwaMBDoBBwGDgPOTOg5LjtkvaXE/UeX1O+It9RGb1f4RcBXwqJm1Ah4EHgohvFZNnQcBLYAJNXw/+wLvVbn/XrJP6kkBLOKhtzyEUFrNY58nj9fHJuDGEMIm4E/J69yVtDBnArPwP+cJIUwLIUwJIZQmre/7gO/V8PrlwC9DCBtCCOs3fzCE8AdgPjAV6AZct5XX2eL7N7PJSWt8vZlV/AJoA6yq8rxVQBv1A9efAlgElgOdttKP2i15vD5WhBDKkq8rAvKLKo+vx0MNM9vTzCYmJ/9W4/3MNQX/shDCNzUc8wegNzC6SlfCFnWy2fcfQhgYQmifPFaRE2uAdlWe1w5YEzSjV70pgEXgLWAD8MOqO82sDfADYFIT1HAv8CGwRwihHXAtUFPLcpvBl9Q/ChgL3GBmHbdyaMX3P7iG95tJ0mJP7Jfsk3pSAEvOCyGswk/CjTazY82smZkVAk8Ci4FHmqCMtsBqYI2Z7Q2ct9njXwC71fE17wKKQwg/Bf4G/L66g0IIX+Hf/z1mdrKZtTWzPDPbH2hd5dCHgcvMrLuZ7QRcDvyxjjVJFQ0auiKSLUIIt5rZCuB2YHc8DJ8FTt/Gn+6pdAUwBrgSmA48ARxZ5fEbgIfMbDv8hFvJtl7MzAYDxwJ9kl2XAe+a2ekhhHGbH598/0uS938YWAsswE/kTU4Ouw//JfBBcv/+ZJ/Uk6n7RkQkDnVBiIhEogAWEYlEASwiEokCWEQkEo2CEAA6deoUCgsLY5chklGmTZu2PITQub7PVwALAIWFhRQXF8cuQySjmNmihjxfXRAiIpEogEVyyapVMHkyLFoEugYgOgWwSLYLAcaPh969oX17OPhgKCyEHXeEq6+Gr7+OXWHOUgCLZLOlS2HQIDjlFMjPhxtvhAkT4J574Igj4JZbYK+94MUXY1eak3QSTiRbLV7s4bt4Mdx9N5xzjodwhfPOg6lTYfhwOPFEePxxGDIkXr05SC1gkWy0ZAkceqi3gP/v/+D8878dvhX694c33oDvfMdbyePHN32tOUwBLJJtNm3yMF2+HCZN8j7fbdl+e++COOggGDYMZsxomjpFASySda680kc6jB0LRUW1e06bNvDUUx7GQ4bA6tWNW6MACmCR7DJxIowaBRdf7K3guujaFZ54Aj76CC64oHHqk29RAItkizVrvK+3d2+47bb6vcZhh8G118Kjj8JLL6W2PtmCAlgkW9xwA3z6Kdx3HzRvXv/XufZa2GMPD/NvalrzUxpCASySDd5/37sehg+HgQMb9lotW/o44fnz4be/TU19Ui0FsEg2uPJKaNcObr45Na931FFw6qnelfHZZ6l5TdmCAlgk0736qg8ju/Za6Li1lefr4aaboLQUfv3r1L2mfIsCWCSTheDzOfToARdemNrX3m03v3ru/vth7tzUvrYACmCRzPbMM/D22/CrX3nfbar9/Of+ur/4RepfWxTAIhkrBPjv//YRC2ee2Tjv0aULXHqpX6I8e3bjvEcOUwCLZKrnn4fp0+Gaa6CgEefVuvRS2G671J3gk39RAItkoorW7847wxlnNO57deoEI0bAuHGwYEHjvleOUQCLZKK//x3eesuHnzVr1vjvd8UVPpvarbc2/nvlEAWwSCa69VZf0eInP2ma9+ve3WdK++MfoaSkad4zByiARTLNrFne/3vhhd4321Quuww2bPCr5CQlFMAimWbUKB8adu65Tfu+e+8Nxx/vAbx+fdO+d5ZSAItkkpISePhhH3bWuXPTv//ll8OyZX5CThpMASySSX7/e+8GGDkyzvsffjjsvz/ccYeWtU8BBbBIpti4Ee69F37wA+8OiMHMw3/2bF/uSBpEASySKZ56yhfZvOiiuHWceqp3f4weHbeOLKAAFskUo0f7ZcfHHBO3jhYt/MKMv/4VFi6MW0uGUwCLZILiYpgyxYee5aXBj+1553kdGpLWIGnwLykiNRo92lcuPuus2JW47t199eT774d162JXk7EUwCLpbvlyX634zDN91Yt0ceGF8NVX8PjjsSvJWApgkXQ3dqwPPTv//NiVfNshh0CfPnD33RqSVk8KYJF0VlbmQ88OPxz23Td2Nd9m5r8Upk/3/mmpMwWwSDp77jlYtAguuCB2JdU74wzvFvnf/41dSUZSAIuks7vvhp12gsGDY1dSvTZtfJa08ePhiy9iV5NxFMAi6WrePF/t+JxzmmbO3/o6/3zYtMlHREidKIBF0tW99/pSQ8OHx65k2/beGwYNgvvu82XspdYUwCLpaN06ePBBH2vbrVvsamp2wQXw6acwcWLsSjKKAlgkHT32mI+xTdeTb5s78UTo2dP7rKXWFMAi6SYED7I+fXysbSYoKPAJ4l9+GT78MHY1GUMBLJJu3nwT3n3XrzQzi11N7f30p9C8uYak1YECWCTdjB4N7dvD6afHrqRudtwRhg6Fhx6C1atjV5MRFMAi6WTJEnj6aTj7bGjdOnY1dXfRRbBmja+eLDVSAIukk9//HsrL02/eh9oqKoL+/b0borw8djVpTwEski6++cbH0p5wAuy2W+xq6u+ii/wikhdeiF1J2lMAi6SLceN8xeFLL41dScP86Ec+dvnOO2NXkvYUwCLpIAQPrL594YgjYlfTMM2b+wiOl1+GDz6IXU1aUwCLpIOXX4aZM33F4UwaerY155wD220Hd90Vu5K0pgAWSQd33gldusBpp8WuJDV22MFX8Hj0USgpiV1N2lIAi8T2wQfw/PN+2XGLFrGrSZ2RI2HjRi1fvw0KYJHYbr3Vx/xmyrwPtbXXXj6P8d13+9hg2YICWCSmRYt8Ucvhw6Fjx9jVpN5VV8HKlfCHP8SuJC0pgEViuuMOP+l22WWxK2kcAwbA977n3+fGjbGrSTsKYJFYSkq8ZXj66T6VY7a66ipYvBgeeSR2JWlHASwSy223+XLz11wTu5LGdeyxcOCBcOONvnSR/IsCWCSGkhK45x748Y/9ZFU2M4MbboCFC9UK3owCWCSG22/3uR9+/vPYlTSN449XK7gaCmCRprZ0qQ/NOu207G/9VqhoBS9YoKkqq1AAizS1X/3KRwTccEPsSprW8cf7qIgbbvBFR0UBLNKk5s71kQ/nnAO9esWupmmZ+UUnn32mOSISCmCRpnTttT5JzfXXx64kjkMP9RWUb74ZVqyIXU10CmCRpvL6677c0BVX+Pppueq3v/VLk3/5y9iVRKcAFmkKpaU+R+4uu8DPfha7mrj23deXXLr3Xl/9OYcpgEWawt13+6xnd94JrVrFria+3/zGp6y84IKcXjtOASzS2D77zPt8jz4aTjopdjXpoX17uOUWmDzZl7HPUQpgkcYUgo942LjRVwrOhtUuUmXYMD8pN3IkLFkSu5ooFMAijenRR2HiRLjpJthjj9jVpJe8PHjgAf/lNGKE/7LKMQpgkcayeDFcfDEcfLDfypZ69fJREc8952GcYxTAIo1h0yYYOtRHPzz4IOTnx64ofV10ka8EfdFFvjBpDlEAizSG666DN9/0q97U9bBteXkwbhy0awc/+lFOLV+kABZJtfHjfa7fc8/1VrDUrFs3eOwx+PBD+MlPcmZomgJYJJUmT4b//E8YONDH/ErtHXmkzxUxfrz/BZEDCmIXIJI15szxVYB79oQJE6Bly9gVZZ7LL4ePPvK5InbeGc47L3ZFjUoBLJIKH37oLTgzP6PfqVPsijKTGYwe7SNIzj/fT16OGBG7qkajLgiRhpoxw8/il5fDa6/ppFtDFRTAU0/5/MHnnOMXsGQpBbBIQzz/vPf3msGrr8I++8SuKDu0aOEzxw0e7MPTRo6EsrLYVaWcAlikPsrK/Oq2E07wiwnefhv+7d9iV5VdKkL40kth1Cg47jhfzimLKIBF6mrBAu9yuO46H7f6+uvQo0fsqrJTfr6PJhkzxj/nPn3g2WdjV5UyCmCR2lq71mc122cfn8f24Yfh8cehTZvYlWW/4cNh2jT/Rfcf/+H9w3Pnxq6qwRTAIjVZtcqHRe26q89j+8MfwqxZPt5Xs5s1nX32galT4fbb4Y03/P6ZZ/q/RYZSAItUp7QUXnkFzjrLr9K65ho48EC/vPixx9TlEEvz5j5WeO5cuOQS7yPed1847DBf7n7lytgV1omFHJwCTrZUVFQUiouLY5cRT1mZX0jx5ps+muHFF+HLL7174cc/9uFQ/frFrlI2t2yZz6I2dizMm+d9xocdBoMG+W2/ftC6daO9vZlNCyEU1fv5CmCBHAjgjRu9dbRsGXzxhQ/0X7QI5s/34J0xA9at82O7doWjjvK+xmOOadQfYEmREOCdd+CZZ+Bvf/Pln8C7iPbc00eo7LEHFBb6lYrduvnCqB07+r9vPbuSFMCSEkWdO4fiwYMb/4229v+tYn/Vx0PYcisv99uyssqttNS3TZtgwwbf1q/3bc0a+Ppr/7o6PXrAXntB797eWurf339g1beb2Vas8Hk5pk/3E6Zz5vgv240btzy2oADatvWtVSvYbjsfAteihXd5NGvmx+TnV2477QR33KEAltQoat48FDfVUulbC7eK/VUfN/v2lpfnt1V/GPLz/YekWbPKH5zttvOtbVvvRmjf3rfOnaFLFw/eHj38WMkN5eX+188nn/htSYn/VbRypf+S/vpr/yto/frKX+QbN/ov9tLSyl/45eW+uvWLLzY4gDUXhLi+fSGbuyBE8vK866Fbt9iV/ItGQYiIRKIAFhGJRH3AAoCZfQ3MiV1HDToBy2MXUYNMqBEyo85MqHGvEELb+j5ZfcBSYU5DTiY0BTMrVo2pkQl1ZkqNDXm+uiBERCJRAIuIRKIAlgpjYhdQC6oxdTKhzqyvUSfhREQiUQtYRCQSBbCISCQK4BxnZsea2Rwzm29mV8eup4KZ9TSzV81slpnNNLNLkv0dzewlM5uX3HZIg1rzzWy6mU1M7u9qZlOTz/QJM2seub72ZvaUmX1oZrPN7KB0+xzNbGTy7zzDzB43s5bp8Dma2QNmVmJmM6rsq/azM/c/Sb3vm1mN85cqgHOYmeUDdwM/APYBTjOzdFnWtxS4PISwDzAAuCCp7WpgUghhD2BScj+2S4DZVe7fAtwZQugFrATOjlJVpbuAF0IIewP74bWmzedoZt2Bi4GiEEJvIB8YSnp8jn8Ejt1s39Y+ux8AeyTbCODeGl89hKAtRzfgIODFKvevAa6JXddWap0AfB+/Wq9bsq8bfgFJzLp6JD+ERwITAcOv3iqo7jOOUN/2wEKSE+5V9qfN5wh0Bz4FOuIXh00EjkmXzxEoBGbU9NkB9wGnVXfc1ja1gHNbxX/8CouTfWnFzAqBA4CpQJcQwufJQ0uBLpHKqjAKuBIoT+7vAHwVQihN7sf+THcFlgEPJt0k95tZa9LocwwhLAFuBz4BPgdWAdNIr8+xqq19dnX+eVIAS1ozszbA08ClIYTVVR8L3syINo7SzE4ASkII02LVUAsFQD/g3hDCAcBaNutuSIPPsQMwGP9lsRPQmi3/7E9LDf3sFMC5bQnQs8r9Hsm+tGBmzfDwHRdC+HOy+wsz65Y83g0oiVUfcDDw72b2MfAnvBviLqC9mVXMsxL7M10MLA4hTE3uP4UHcjp9jkcBC0MIy0IIm4A/459tOn2OVW3ts6vzz5MCOLe9A+yRnG1ujp/4+EvkmgA/owyMBWaHEO6o8tBfgGHJ18PwvuEoQgjXhBB6hBAK8c/ulRDC6cCrwMnJYbFrXAp8amZ7JbsGAbNIo88R73oYYGatkn/3ihrT5nPczNY+u78AZyajIQYAq6p0VVQvVse7tvTYgOOAucBHwHWx66lS1yH4n3bvA+8m23F4H+skYB7wMtAxdq1JvYcDE5OvdwPeBuYD44EWkWvbHyhOPstngQ7p9jkCvwI+BGYAjwAt0uFzBB7H+6U34X9NnL21zw4/AXt38rP0AT6qY5uv36BLkc3sWPxPrnzg/hDCzfV+MRGRHFPvAE7GkM7FhwYtxv+cPS2EMCt15YmIZK+GTMj+XWB+CGEBgJn9CT+TudUA7tSpUygsLGzAW0pDffWVL/7as+e390+bNm15CKFzxf3v5/1IszSJVPFS+fitLOddfw0J4OrGvPXf/CAzG4FfFcLOO+9MsVbejerKK2H06C0XQDazRXEqEsldjT4KIoQwJoRQFEIo6ty5c81PkEa1aRM0axa7ChGBhgVwWo8hlept3AjNo04NIyIVGhLAaTuGVLZu3Tpo3Tp2FSICDegDDiGUmtmFwIv4MLQHQggzU1aZNIrVq6FNm9hViAg0cFn6EMJzwHMpqkWawIoVsMMOsasQEdClyDlnyRLo2jV2FSICCuCcUlYGH38Mu+8eu5KIzHwTSQMK4BwyYwaUlkLv3rErERFoYB+wZJa//91vDzssbh1RVVx6X9EKtqQNEsq//bhIE1ALOIc8+ijss8+WlyGLSBxqAeeIKVPgnXf8MmShsqUbyvw2Lx8Aa+4/EmFTshJOeVlTVyY5RC3gHLBxI5x3HnTpAmeeGbsaEamgFnCWCwGuuw7efReefRbatYtdUZpKWrphg99agf9o2HZ+2WDYuMlvN22MUJxkK7WAs1gIPvvZ7bd7C3jw4NgViUhVagFnqVWr4OKL4eGH4YIL4H/+J3ZFmSWUln7rNq9lS7/d3i8jDOvWA1C+bl2E6iRbqAWchZ591kc7PPooXH+9n3jL07+0SNpRCziLvPUW/Pa38Ne/Qt++HsTf+U7sqrJD+Tff+BfJbX7SmV6wi4/pC6vXAFC2cmXTFycZS+2iDFdWBs88AwcfDAMHwj/+ATfd5CteKHxF0ptawBlq9mx44gkYNw7mz4fCQu/n/a//0nSTTaFs9Wr/IrnNT1Z7yeu7t98u+wqA0s+XNn1xkjEUwBlk7lwP3Sef9HkdzPyy4htvhB/+EAr0rymSUfQjm8bWrIE33oCXX4aXXoIPPvDQPeQQP7E2ZAh06xa7SgEoW7bMv0hu8wp39v1H9AOg+WJvEZfNW9D0xUnaUgCnkU2bYOpUD9xJk/zy4dJSX8Pt4INh1Cg4+WTo3j12pSKSCgrgiEpKPHArtrfegrVrvZV74IFw+eVw1FEevtttF7taqYvSjz8BID+5Zd+9APh66AAA2s37GoAwTat45TIFcBPZsMEvB54yxcN2yhRYuNAfy8+HPn1g2DAYNAgOPxw6doxarog0AQVwI1izBt5/3wN3+nS/ff99nxQHvAthwAC/PLh/f2/taqXi7FY2cw4AbZMGb/nA/QD4YuRAADq9twGAglemNX1xEo0CuIG++OLbQTt9OsybVznbYceOcMABcMklHrb9+0OPHnFrFpH0oACupfXrYdYsH4lQdVtaZZhnYaGH7emnw/77+9c9emgJMtmSTX4PgK6T/f6G4/yqmUW3HQRAt7d8hY5Wf57a9MVJk1EAb6asDBYs2DJo58+H8mTVmpYtfa6FY46pDNr99oP27ePWLiKZJacDuKRky6CdORMqJrgy8xWE+/SBoUP9tk8f6NXLT5yJpEqL594BYPfn/P7qH/toiYV/6gvA9v/nJwk6PvBW0xcnjSYnAnjdOg/WzcO2pKTymM6dPVyHD68M2n331ckxEWk8WRfAX37pJ8L++c/K27lzK0+KtWzpy7Iff3xl0Pbp48v1iKSLdo9NSW79fsn5Plqi9es+58RHz+wBQNc7Jzd9cZIyGRvAIcBnn20Ztp98UnlMz57Qr9+3uw92313dByKSHjImgMvLvRvh9dd9foTXX4fPP698fM894aCDfPWHAw7wrVOnePWKpNKO93hLd+09fr/sam8BnzDT5x++90/HA9DzN2oRZ5K0DeBNm2DatMqwffNNqJjrunt3v1pswABv4e63H7RtG7VcEZE6S6sADsEnFH/4YZ9ysWLK1T339OkWDz3Up18sLNTYWslt3W/2lu7EmzsAUPYbP8lxy0IfN3zKE5cCsOvVGjWRzmoMYDPrCTwMdAECMCaEcJeZdQSeAAqBj4FTQgj1Wo9lwQJ45BEP3gULfOTBkCFw4ok+9WLXrvV5VRGR9FabFnApcHkI4Z9m1haYZmYvAWcBk0IIN5vZ1cDVwFV1LeB3v4MrrvAW7aBBcMMN3trV8C+R2iv8hbd0r/pFfwDKb/P9Ty/20RT7Jy3i3S+f0vTFyVbVGMAhhM+Bz5Ovvzaz2UB3YDBweHLYQ8Br1DGAx4zx8B0yBO6800ctiIjkijr1AZtZIXAAMBXokoQzwFK8i6LWFi6Ec8+FoiJ47DGfdFxEUmP3n3mLeMjP/Iq6cKfvf/Gzd/3xP50LQK/L1CKOqdarIptZG+Bp4NIQwuqqj4UQAt4/XN3zRphZsZkVL6tYtgWfpOaII3z87nPP1a94EZFMVqsANrNmePiOCyH8Odn9hZl1Sx7vBpRU99wQwpgQQlEIoahzsnIsQLNm8Oyz3gIeMgSOO84XnPzmmwZ9PyJSjV4jp9Br5BSO2Wl/jtlpfyyABe8jfnrxFD667SA+SmZik6ZTYwCbmQFjgdkhhDuqPPQXYFjy9TBgQl3fvG1beOEFuOoqn5th6FAf8TBihI/7rZh9TEQkG1kI1fYcVB5gdgjwBvABUBGJ1+L9wE8COwOL8GFoX27rtYqKikJxcXG1j5WVwWuvwUMPwdNP+wQ6HTr4MLTDDvMxwP36ectZUs/MpoUQiirufz/vR9v+jyFZZeHN3vp98tRRAJz6qI+aqBhdIfBS+fiUX31Qm1EQ/wC29saDUlVIfr4PQxs0CO65ByZMgFdf9Svh/vpXP6ZVK7/cuOKCjO9+V8PVRCRz1dgCTqVttYC3ZelSv0KuYh6I997zq+bMYK+9fN6Hfv0q54DQgpZ1pxawVPXJ9T772vmn/g2A+x7xuSYqrsDLRVFawOmga1c4+WTfAL76CiZPhnfe8VEU//gHPP545fG77FIZyP36+aoVO+2ky5dFJL1kRABvrn17HzVx3HGV+5Yv9zCu2P75Tx9lUdHA79ChckrKvn39tndvTeIjUp2df53MNfFrn2sif6Tvr5iP+OPHewHQ+V71ETdERnRB1NfXX3t3xXvv+SiL99+HGTN8f4XCwm9PzN6nj0/+k2sn+9QFIXXx5U/8pN2qo9cCsMNfWgGVE8lno5ztgqivtm19FMUhh1TuCwEWLdpyeaLnn4fSUj+meXPYe+8tg1krHItIKmV1C7guNmyAOXO8lVw1mBcvrjymfXvvtti8K6Ndu3h1p4pawNIQa4f4JEBLB/ilBT1e9dZMxWKj2UAt4EbUooUHat++396/cqV3W1QN5XHjKucqBl8luWJ5+orbrl3VWhaRbVMLuB5CgE8/9dbyu+/6Nn26z2VcYccdKwN5//3hwAM9qNM1lNUCllQqPfJAAJbv14IuU72f2Ca/F7OkBlMLOE2Ywc47+3bCCZX7V63yUJ4+vTKU77jDl1cCH5/cv79vAwb4hSQdOsT5HkQkPgVwCm2/vV+ld+ihlfs2boRZs6C4GKZOhSlTfP6Lij889tzTw7h/fxg40LtA8mo9R51Ieip4ZRoAXV8BO3BfAFYP9akx23/gC+eUzZwTp7g0ogBuZM2bV3ZD/PSnvm/1ag/kKVM8lF94wZdjAl/J+Ygj4Kij/LLs3XZL324LEWkYBXAE7drBkUf6BpVD415/HSZN8m38eH+ssLByjoxjj1WXhWSeMG0mAG2nJTv22A2AsiP6AdBigc8TXrro0yavLTYFcBow86AtLIQzz/RAnjOnMoyffhrGjvWLQ44+Gk45BQYP9i4PEclcCuA0ZOYXguy9N1xwgU/V+c47HsRPPgl/+5t3bRx7rIfxkCHQsmXsqkVqp2yeDxfKn5fs6L4TAHl99wbAPl/hx1VZQSdb6XRPBsjP9xN1t90GH38Mb73lwTxtGpxxhk8+9JvfwIoVsSsVkbpQAGcYMw/jO+6ATz6Bl1/2McbXX++rSl94ofcni2SK0iWfUbrkM8rf/5Dy9z/0OQFKSynYpScFu/Qkr21b8rJ01iwFcAbLy/OTc889V7mk05gxsM8+MGqUd12ISPpSAGeJ3r3hgQdg3jw4/HAYOdLHFc+cGbsykbopW7mSspUrKV30qY+MKC+H8nLyO+1AfqcdyGvZkrwsOemhAM4yu+wCEyfCY4/5pdEDB/oCpyKSfhTAWcgMTjvNJ6Xv2tWHrk2aFLsqkfopX7uW8rVrKVu+grLlKwghEEIgr3Vr8lq3xgoKsILMHNClAM5iPXv6xR277ebLOS1ZErsiEalKAZz3yuHQAAAGyklEQVTlunSBZ57x+Y6HD6+cg0IkU4UNGwgbNvyrZVzBWrTAWrTwPwEz5Pp9BXAO6NULbrrJV/149dXY1YhIBQVwjjj3XJ8O8557YlciklqhtNS3pGX8L3n5vqUxBXCOaNkShg2DCROgyl9tIhKRAjiHHH20X2Q0JXsXrhXxEx0hQHmZbxV9wmnYN6wAziEH+UrivP123DpExGXm4Dmpl+2391ERH30UuxKRJpTGQ3/UAs4xhYU+iY+IxKcAzjGdOmnaSpF0oQDOMR06wMqVsasQEahDAJtZvplNN7OJyf1dzWyqmc03syfMrHnjlSmp0qoVrF8fuwoRgbq1gC8BZle5fwtwZwihF7ASODuVhUnjaNlSASySLmoVwGbWAzgeuD+5b8CRwFPJIQ8BJzVGgZJazZrBpk2xqxARqH0LeBRwJVCe3N8B+CqEUJrcXwx0r+6JZjbCzIrNrHhZDiyyl+4KCvxiDBGJr8YANrMTgJIQwrT6vEEIYUwIoSiEUNS5c+f6vISkUF6eLzAgIvHV5kKMg4F/N7PjgJZAO+AuoL2ZFSSt4B6AZpvNAApgkfRRYws4hHBNCKFHCKEQGAq8EkI4HXgVODk5bBgwodGqlJQxS+sLg0RySkPGAV8FXGZm8/E+4bGpKUkakwJYJH3UaS6IEMJrwGvJ1wuA76a+JGlMaTYZlEhO05VwIiKRKIBzjFrAIulDASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikdQqgM2svZk9ZWYfmtlsMzvIzDqa2UtmNi+57dDYxYqIZJPatoDvAl4IIewN7AfMBq4GJoUQ9gAmJfdFRKSWagxgM9seOAwYCxBC2BhC+AoYDDyUHPYQcFJjFSkiko1q0wLeFVgGPGhm083sfjNrDXQJIXyeHLMU6FLdk81shJkVm1nxsmXLUlO1iEgWqE0AFwD9gHtDCAcAa9msuyGEEIBQ3ZNDCGNCCEUhhKLOnTs3tF4RkaxRmwBeDCwOIUxN7j+FB/IXZtYNILktaZwSRUSyU40BHEJYCnxqZnsluwYBs4C/AMOSfcOACY1SoYhIliqo5XEXAePMrDmwAPgvPLyfNLOzgUXAKY1ToohIdqpVAIcQ3gWKqnloUGrLERHJHboSTkQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiERSqwA2s5FmNtPMZpjZ42bW0sx2NbOpZjbfzJ4ws+aNXayISDapMYDNrDtwMVAUQugN5ANDgVuAO0MIvYCVwNmNWaiISLapbRdEAbCdmRUArYDPgSOBp5LHHwJOSn15IiLZq8YADiEsAW4HPsGDdxUwDfgqhFCaHLYY6F7d881shJkVm1nxsmXLUlO1iEgWqE0XRAdgMLArsBPQGji2tm8QQhgTQigKIRR17ty53oWKiGSb2nRBHAUsDCEsCyFsAv4MHAy0T7okAHoASxqpRhGRrFSbAP4EGGBmrczMgEHALOBV4OTkmGHAhMYpUUQkO9WmD3gqfrLtn8AHyXPGAFcBl5nZfGAHYGwj1ikiknUKaj4EQgi/BH652e4FwHdTXpGISI7QlXAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwDlmwAC45JLYVYgIgIUQmu7NzJYBi5rsDaUudgkhdI5dhEguadIAFhGRSuqCEBGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKR/D+iA6AGU5wC6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% EMD\n",
    "\n",
    "G0 = ot.emd(a, b, M)\n",
    "\n",
    "pl.figure(3, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, G0, 'OT matrix G0')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve Sinkhorn\n",
    "--------------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It.  |Err         \n",
      "-------------------\n",
      "    0|8.187970e-02|\n",
      "   10|3.460174e-02|\n",
      "   20|6.633335e-03|\n",
      "   30|9.797798e-04|\n",
      "   40|1.389606e-04|\n",
      "   50|1.959016e-05|\n",
      "   60|2.759079e-06|\n",
      "   70|3.885166e-07|\n",
      "   80|5.470605e-08|\n",
      "   90|7.702918e-09|\n",
      "  100|1.084609e-09|\n",
      "  110|1.527180e-10|\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcVNWZ//HP0xv7KoRdW0RhFFyQBFxHxTWaIRON0ZgfmBjc15jRqJlEJxMTjeMyRE2MxGBERsG4hLhEjUajgoIYRRBEFAVBmn2nt/P747lFN01DN72dWr7v1+u+iqq6feuhoL/99Klzz7UQAiIi0vLyYhcgIpKrFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWCQyMzvKzOY1w3HPMbO/1nPfc83sH7v7nDSOAlgyXhIQ75nZJjNbZmb3mlnn5LnfmNmGZCs1s7Jq959pgdqCmQ3Y1T4hhFdDCAMbePwjzex1M1trZqvM7DUz+3Jy3IkhhBMbclxpGQpgyWhmdjVwC/AfQCdgBLAX8LyZFYUQLgwhtA8htAduBh5J3Q8hnBKvcmdmBY342o7AVGAc0BXoA9wEbG2a6ppeY/6+2UgBLBkrCaCbgMtCCM+GEMpCCJ8AZwLFwHcacMxjzGyxmV1jZsvNbKmZfd3Mvmpm85Mu8/pq+3/FzN4wszXJvr82s6LkuVeS3f6ZdNzfqnb8a81sGfBA6rHka/ZJXmNocr+3mZWY2TG1lLsfQAhhUgihIoSwOYTw1xDCu8nXbjd0kHTjF5rZh0m9d5uZ7eR9+JWZ/cPMOlV77DYzW21mH5vZKdUe721mTyV1LzCzsdWeu9HMppjZQ2a2Djg3eexRM3vQzNab2ftmNmw3/6myggJYMtnhQGvgT9UfDCFsAJ4GTmjgcXsmx+0D/AT4HR7mhwJHAf9pZnsn+1YAVwHdgMOAkcDFSR1HJ/sclHTcj1Q7fle8Uz+/Ru0fAdcCD5lZW+ABYEII4eVa6pwPVJjZBDM7xcy61OPvdhrwZeBA/AfVSdWfNLM8M/td8vyJIYS1yVPDgXnJ3/NWYHy18P4/YDHQGzgDuNnMjqt22FHAFKAzMDF57N+Sr+sMPAX8uh61Zx0FsGSybsCKEEJ5Lc8tTZ5viDLg5yGEMjwkugF3hRDWhxDeB+YABwGEEGaGEKaFEMqT7vu3wL/WcfxK4KchhK0hhM01nwwh/A5YAEwHegE31HaQEMI64Egg4D8kSpJOtMcuXvuXIYQ1IYRPgZeAg6s9VwhMwn84fC2EsKnac4tCCL8LIVQAE5K6ephZP+AI4NoQwpYQwjvA/cDoal/7RgjhiRBCZbW/7z9CCE8nx/sjyfuZaxTAkslWAN12Mq7YK3m+IVYmwQCQCowvqj2/GWgPYGb7mdnU5MO/dfg4c13BXxJC2FLHPr8DBgPjQgg7HdMNIcwNIZwbQuib7N8buHMXx11W7c+bUn+PxAC8W70phFC6s6+rFsztk9dbFUJYX23fRfhvDymf1aOO1rk4PqwAlkz2Bv6B0zeqP2hm7YFTgBdboIZ7gQ+AfUMIHYHrgVrHVavZ5RKESf13AuOBG82sa30KCSF8APwBD+KGmAt8F3jGzOo7K+NzoKuZdaj22J7AkuqlNbCerKcAloyVjE/eBIwzs5PNrNDMioFH8THJP7ZAGR2AdcAGMxsEXFTj+S+A/rt5zLuAGSGE7wN/AX5T205mNsjMrjazvsn9fsDZwLTdfL1tQgiT8B8iL5jZPvXY/zPgdeAXZtbazA4EzgMeamgNuUQBLBkthHArHhi34UE4Hf+Vd+SufnVvQj8Evg2sx4cNHqnx/I3AhGTWwZl1HczMRgEnUxXkPwCGmtk5tey+Hv9wbLqZbcSDdzZwdQP+HtuEECYA/wX8LfmBVpez8VknnwOP4+PbLzSmhlxhWpBdRCQOdcAiIpEogEVEIlEAi4hEogAWEYkk5yY+S+26desWiouLY5chklFmzpy5IoTQvaFfrwAWAIqLi5kxY0bsMkQyipktaszXawhCRCQSBbBILlm7Fl5/HRYtAp0DEJ0CWCTbhQCTJ8PgwdC5MxxxBBQXw5e+BD/6EaxfX+chpHkogEWy2bJlMHIknHkm5OfDz38OTz4J99wDxx4Lt9wCAwfCc8/FrjQn6UM4kWy1eLGH7+LFcPfdcMEFHsIpF10E06fD2LHwta/BpElw+unx6s1B6oBFstGSJXDUUd4B//WvcPHF24dvyvDh8Oqr8OUve5c8eXLL15rDFMAi2aaszMN0xQp48UUf892VTp18COKww2DMGJg9u2XqFAWwSNa55hqf6TB+PAyr57Uu27eHKVM8jE8/Hdata94aBVAAi2SXqVPhzjvh8su9C94dPXvCI4/ARx/BJZc0T32yHQWwSLbYsMHHegcPhl/9qmHHOPpouP56eOgheP75pq1PdqAAFskWN94In30Gv/0tFBU1/DjXXw/77uthvqWua4dKYyiARbLBu+/60MPYsXD44Y07VuvWPk94wQL4xS+apj6plQJYJBtccw107Ai//GXTHO/44+Fb3/KhjM8/b5pjyg4UwCKZ7qWXfBrZ9ddD13pdwb5+br4Zysvhv/6r6Y4p21EAi2SyEHw9h7594dJLm/bY/fv72XP33w/z5zftsQVQAItktscfhzffhJtu8rHbpvbjH/tx//M/m/7YogAWyVghwH//t89YGD26eV6jRw+48ko/RXnu3OZ5jRymABbJVM88A7NmwXXXQUEzrqt15ZXQpk3TfcAn2yiARTJRqvvdc0/4znea97W6dYPzz4eJE2HhwuZ9rRyjABbJRH//O7zxhk8/Kyxs/tf74Q99NbVbb23+18ohCmCRTHTrrX5Fi+99r2Ver08fXyntD3+A5ctb5jVzgAJYJNPMmePjv5de6mOzLeUHP4CtW/0sOWkSCmCRTHPnnT417MILW/Z1Bw2CU0/1AN68uWVfO0spgEUyyfLl8OCDPu2se/eWf/2rr4aSEv9AThpNASySSX7zGx8GuOqqOK9/zDFw8MFw++26rH0TUACLZIrSUrj3XjjlFB8OiMHMw3/uXL/ckTSKAlgkU0yZ4hfZvOyyuHV861s+/DFuXNw6soACWCRTjBvnpx2fdFLcOlq18hMz/vxn+PjjuLVkOAWwSCaYMQOmTfOpZ3lp8G170UVeh6akNUoa/EuKSJ3GjfMrF597buxKXJ8+fvXk+++HTZtiV5OxFMAi6W7FCr9a8ejRftWLdHHppbBmDUyaFLuSjKUAFkl348f71LOLL45dyfaOPBKGDIG779aUtAZSAIuks4oKn3p2zDFwwAGxq9memf9QmDXLx6dltymARdLZ00/DokVwySWxK6ndd77jwyK//nXsSjKSAlgknd19N/TuDaNGxa6kdu3b+yppkyfDF1/EribjKIBF0tWHH/rVji+4oGXW/G2oiy+GsjKfESG7RQEskq7uvdcvNTR2bOxKdm3QIBg5En77W7+MvdSbAlgkHW3aBA884HNte/WKXU3dLrkEPvsMpk6NXUlGUQCLpKOHH/Y5tun64VtNX/sa9OvnY9ZSbwpgkXQTggfZkCE+1zYTFBT4AvEvvAAffBC7moyhABZJN6+9Bu+842eamcWupv6+/30oKtKUtN2gABZJN+PGQefOcM45sSvZPV/6Epx1FkyYAOvWxa4mIyiARdLJkiXw2GNw3nnQrl3sanbfZZfBhg1+9WSpkwJYJJ385jdQWZl+6z7U17BhMHy4D0NUVsauJu0pgEXSxZYtPpf2tNOgf//Y1TTcZZf5SSTPPhu7krSnABZJFxMn+hWHr7wydiWN881v+tzlO+6IXUnaUwCLpIMQPLAOPBCOPTZ2NY1TVOQzOF54Ad57L3Y1aU0BLJIOXngB3n/frzicSVPPduaCC6BNG7jrrtiVpDUFsEg6uOMO6NEDzj47diVNY489/AoeDz0Ey5fHriZtKYBFYnvvPXjmGT/tuFWr2NU0nauugtJSXb5+FxTAIrHdeqvP+c2UdR/qa+BAX8f47rt9brDsQAEsEtOiRX5Ry7FjoWvX2NU0vWuvhdWr4Xe/i11JWlIAi8R0++3+odsPfhC7kuYxYgT867/637O0NHY1aUcBLBLL8uXeGZ5zji/lmK2uvRYWL4Y//jF2JWlHASwSy69+5Zebv+662JU0r5NPhkMPhZ//3C9dJNsogEViWL4c7rkHvv1t/7Aqm5nBjTfCxx+rC65BASwSw223+doPP/5x7EpaxqmnqguuhQJYpKUtW+ZTs84+O/u735RUF7xwoZaqrEYBLNLSbrrJZwTceGPsSlrWqaf6rIgbb/SLjooCWKRFzZ/vMx8uuAAGDIhdTcsy85NOPv9ca0QkFMAiLen6632Rmp/8JHYlcRx1lF9B+Ze/hJUrY1cTnQJYpKW88opfbuiHP/Trp+WqX/zCT03+6U9jVxKdAlikJZSX+xq5e+0F//EfsauJ64AD/JJL997rV3/OYQpgkZZw992+6tkdd0DbtrGrie9nP/MlKy+5JKevHacAFmlun3/uY74nnghf/3rsatJD585wyy3w+ut+GfscpQAWaU4h+IyH0lK/UnA2XO2iqYwZ4x/KXXUVLFkSu5ooFMAizemhh2DqVLj5Zth339jVpJe8PPj97/2H0/nn+w+rHKMAFmkuixfD5ZfDEUf4rexowACfFfH00x7GOUYBLNIcysrgrLN89sMDD0B+fuyK0tdll/mVoC+7zC9MmkMUwCLN4YYb4LXX/Kw3DT3sWl4eTJwIHTvCN7+ZU5cvUgCLNLXJk32t3wsv9C5Y6tarFzz8MHzwAXzvezkzNU0BLNKUXn8d/t//g8MP9zm/Un/HHedrRUye7L9B5ICC2AWIZI158/wqwP36wZNPQuvWsSvKPFdfDR995GtF7LknXHRR7IqalQJYpCl88IF3cGb+iX63brErykxmMG6czyC5+GL/8PL882NX1Ww0BCHSWLNn+6f4lZXw8sv60K2xCgpgyhRfP/iCC/wEliylABZpjGee8fFeM3jpJdh//9gVZYdWrXzluFGjfHraVVdBRUXsqpqcAlikISoq/Oy2007zkwnefBP+5V9iV5VdUiF85ZVw553w1a/65ZyyiAJYZHctXOhDDjfc4PNWX3kF+vaNXVV2ys/32ST33efv85Ah8MQTsatqMgpgkfrauNFXNdt/f1/H9sEHYdIkaN8+dmXZb+xYmDnTf9D9+7/7+PD8+bGrajQFsEhd1q71aVF77+3r2H7jGzBnjs/31epmLWf//WH6dLjtNnj1Vb8/erT/W2QoBbBIbcrL4W9/g3PP9bO0rrsODj3UTy9++GENOcRSVORzhefPhyuu8DHiAw6Ao4/2y92vXh27wt1iIQeXgJMdDRs2LMyYMSN2GfFUVPiJFK+95rMZnnsOVq3y4YVvf9unQw0dGrtKqamkxFdRGz8ePvzQx4yPPhpGjvTboUOhXbtme3kzmxlCGNbgr1cAC+RAAJeWendUUgJffOET/RctggULPHhnz4ZNm3zfnj3h+ON9rPGkk5r1G1iaSAjw1lvw+OPwl7/45Z/Ah4j2289nqOy7LxQX+5mKvXr5hVG7dvV/3wYOJSmApUkM6949zBg1qvlfaGf/31KPV38+hB23ykq/raio2srLfSsrg61bfdu82bcNG2D9ev9zbfr2hYEDYfBg75aGD/dvWI3tZraVK31djlmz/APTefP8h21p6Y77FhRAhw6+tW0Lbdr4FLhWrXzIo7DQ98nPr9p694bbb1cAS9MYVlQUZrTUpdJ3Fm6px6s/b7b9lpfnt9W/GfLz/ZuksLDqG6dNG986dPBhhM6dfeveHXr08ODt29f3ldxQWem//Xz6qd8uX+6/Fa1e7T+k16/334I2b676QV5a6j/Yy8urfuBXVvrVrZ97rtEBrLUgxB14IGTzEIRIXp4PPfTqFbuSbTQLQkQkEgWwiEgkGgMWAMxsPTAvdh116AasiF1EHTKhRsiMOjOhxoEhhA4N/WKNAUvKvMZ8mNASzGyGamwamVBnptTYmK/XEISISCQKYBGRSBTAknJf7ALqQTU2nUyoM+tr1IdwIiKRqAMWEYlEASwiEokCOMeZ2clmNs/MFpjZj2LXk2Jm/czsJTObY2bvm9kVyeNdzex5M/swue2SBrXmm9ksM5ua3N/bzKYn7+kjZlYUub7OZjbFzD4ws7lmdli6vY9mdlXy7zzbzCaZWet0eB/N7PdmttzMZld7rNb3ztz/JvW+a2Z1rl+qAM5hZpYP3A2cAuwPnG1m6XJZ33Lg6hDC/sAI4JKkth8BL4YQ9gVeTO7HdgUwt9r9W4A7QggDgNXAeVGqqnIX8GwIYRBwEF5r2ryPZtYHuBwYFkIYDOQDZ5Ee7+MfgJNrPLaz9+4UYN9kOx+4t86jhxC05egGHAY8V+3+dcB1sevaSa1PAifgZ+v1Sh7rhZ9AErOuvsk34XHAVMDws7cKanuPI9TXCfiY5AP3ao+nzfsI9AE+A7riJ4dNBU5Kl/cRKAZm1/XeAb8Fzq5tv51t6oBzW+o/fsri5LG0YmbFwCHAdKBHCGFp8tQyoEekslLuBK4BKpP7ewBrQgjlyf3Y7+neQAnwQDJMcr+ZtSON3scQwhLgNuBTYCmwFphJer2P1e3svdvt7ycFsKQ1M2sPPAZcGUJYV/254G1GtHmUZnYasDyEMDNWDfVQAAwF7g0hHAJspMZwQxq8j12AUfgPi95AO3b8tT8tNfa9UwDntiVAv2r3+yaPpQUzK8TDd2II4U/Jw1+YWa/k+V7A8lj1AUcA/2ZmnwD/hw9D3AV0NrPUOiux39PFwOIQwvTk/hQ8kNPpfTwe+DiEUBJCKAP+hL+36fQ+Vrez9263v58UwLntLWDf5NPmIvyDj6ci1wT4J8rAeGBuCOH2ak89BYxJ/jwGHxuOIoRwXQihbwihGH/v/hZCOAd4CTgj2S12jcuAz8xsYPLQSGAOafQ+4kMPI8ysbfLvnqoxbd7HGnb23j0FjE5mQ4wA1lYbqqhdrIF3bemxAV8F5gMfATfErqdaXUfiv9q9C7yTbF/Fx1hfBD4EXgC6xq41qfcYYGry5/7Am8ACYDLQKnJtBwMzkvfyCaBLur2PwE3AB8Bs4I9Aq3R4H4FJ+Lh0Gf7bxHk7e+/wD2DvTr6X3sNndezy+I06FdnMTsZ/5coH7g8h/LLBBxMRyTENDuBkDul8fGrQYvzX2bNDCHOarjwRkezVmAXZvwIsCCEsBDCz/8M/ydxpAHfr1i0UFxc34iWlsdas8Yu/9uu3/eMzZ85cEULonrp/Qt43tUqTSDXPV07eyeW8G64xAVzbnLfhNXcys/Pxs0LYc889maEr70Z1zTUwbtyOF0A2s0VxKhLJXc0+CyKEcF8IYVgIYVj37t3r/gJpVmVlUFgYuwoRgcYFcFrPIZXalZZCUdSlYUQkpTEBnLZzSGXnNm2Cdu1iVyEi0Igx4BBCuZldCjyHT0P7fQjh/SarTJrFunXQvn3sKkQEGnlZ+hDC08DTTVSLtICVK2GPPWJXISKgU5FzzpIl0LNn7CpEBBTAOaWiAj75BPbZJ3YlEZjt3ibSAhTAOWT2bCgvh8GDY1ciItDIMWDJLH//u98efXTcOtKC1dF7pJrgULn9441YO0WkJnXAOeShh2D//Xc8DVlE4lAHnCOmTYO33vLTkHNSqnNNje+mOtukE7Y82+4+qfuVNTre5OtC6vHUcdQZSwOoA84BpaVw0UXQoweMHh27GhFJUQec5UKAG26Ad96BJ56Ajh1jVxRZzU54Jyz1fFHt3yKWHCdUVPgDlTXuqzOWelAHnMVC8NXPbrvNO+BRo2JXJCLVqQPOUmvXwuWXw4MPwiWXwP/+b+yK0kyqMw0VyV3veC2f7e8nu1th8q2Sn+yQdMipTpik8w0VSedbVrbdfXXGUht1wFnoiSd8tsNDD8FPfuIfvOXpX1ok7agDziJvvAG/+AX8+c9w4IEexF/+cuyqMkSNMd1t04BTzyc/wSzpgC21qHJB8i2UzJqw1OyI8nL/+lQnvLXUH0/dTz2/rTNWR5yL1BdluIoKePxxOOIIOPxw+Mc/4Oab/YoXCl+R9KYOOEPNnQuPPAITJ8KCBVBc7OO83/2ulptslJ3MbthBXtIJJx1waJWscp+MFYf8pCOu8K+3Uu98bYt3wmHz5uR2i9+WpsaMNVacSxTAGWT+fA/dRx/1dR3M/LTin/8cvvGNqt+GRSQz6Fs2jW3YAK++Ci+8AM8/D++956F75JH+wdrpp0OvXrGrzFI1Z0mUJh1pqiOuTHWofps6ky4kHXBlW++IK1t5pxxSY8Tl/vX5m73jzVvvHXDexqQj3rjJb1MdssaKs5oCOI2UlcH06R64L77opw+Xl/s13I44Au68E844A/r0iV2piDQFBXBEy5d74Ka2N96AjRu9yz30ULj6ajj+eA/fNm1iV5vjUmPD5duP1ealOtNkvm9qFkReMk+4stD/4co6+LdaeZukU85rBUB+mV+gr3CDH6do9Vb/+jUb/XjrNvj+G/z+trHipA51xJlNAdxCtm7104GnTfOwnTYNPv7Yn8vPhyFDYMwYGDkSjjkGunaNWq6ItAAFcDPYsAHefdcDd9Ysv333XV8UB3wIYcQIPz14+HDvdnWl4gxRY2y4cmvS+aY64OQfOS+Z71tYVp7s71NTKou8893c0Tvh8nY+EzQksyryt/j84tZr/D9Em5IuABSt8LHh/JXrfP8N3hlXpmZRpF6nsqIJ/pLSUhTAjfTFF9sH7axZ8OGHVd+nXbvCIYfAFVd42A4fDn37xq1ZRNKDArieNm+GOXN8JkL1bdmyqn2Kiz1szzkHDj7Y/9y3ry4xltVSY8Nlyfze1JoQyeyFvOQMuKLkNn9rJwCswseGNxR657ulu/8n2ZTMalmf/Kcp2Ogdc+sVfttumS9n1/Zz73wLv1jrr7c26YyTWRSVpakxYs0nTmcK4BoqKmDhwh2DdsGCqplHrVv7WgsnnVQVtAcdBJ07x61dRDJLTgfw8uU7Bu3778MmbyIw8ysIDxkCZ53lt0OGwIABVYtiiWyncvt5w5WptSW2JrMbkjHb9hv9p3XBxg4ArC/1ecMb+nnnu6VXMqbb1zvn9QP97pp1vl/rz32MuMOnbZPbPQBotdQ74fwVa7yOjanZE9t36OqI00NOBPCmTR6sNcN2+fKqfbp393AdO7YqaA84QB+OiUjzyboAXrXKPwh7++2q2/nzq37gt27tl2U/9dSqoB0yxC/XI9Jkts0b3v5MNktmK6Q64TbrvEMtWuOzHVqt9Z/4azf5t+aGffxwXft4R9uvz1IAygf57IlPVvl8xZWLvJPusLA7AJ0Xeofd5rP1AOQvXw1A5fpkXnHSkasjjitjAzgE+PzzHcP200+r9unXD4YO3X74YJ99NHwgIukhYwK4stKHEV55xddHeOUVWLq06vn99oPDDvOrPxxyiG/dusWrV2Q7NWZLVNQcG05mL3Ra5bMk2iz3jnhtSWsA1gz0Md4FA32e8BF9FwJwfPe5AGwa4LMkXj+kPwDvL/Tz1dvN8w65y3w/brtPvAPe1hGvSWZRqCOOIm0DuKwMZs6sCtvXXoPV/n+GPn38bLERI7zDPegg6NAharkiIrstrQI4BF9Q/MEHfcnFdf6BLvvt58stHnWUL79YXKy5tZLhUrMltvptxbaxYV8FrXCdj912W+EdbLul3sGu/txv/3rAEACWDPKx3jN7vgXATXu9B8DKvj6WPPXAgwF49sP9ASh437++yzw/M6/DQu+085epI46hzgA2s37Ag0AP/Aot94UQ7jKzrsAjQDHwCXBmCGF1Q4pYuBD++EcP3oULfebB6afD177mSy/27NmQo4qIpDcLdfxkM7NeQK8Qwttm1gGYCXwdOBdYFUL4pZn9COgSQrh2V8caNmxYmDFjxnaP/c//wA9/6B3tyJEwerR3u5r+1bLMbGYIYVjq/gl531TLE0PqVzvzWQ55rX1s1zp4xxp6eEe8cW8/I27VIO+hNg/xzvnUQbMBGLPHawAMLPT5yJ8nneykNX6dqkc+HOrHe9eP0/UD36/DR8msCXXEO3i+cnKT/95dZwccQlgKLE3+vN7M5gJ9gFHAMcluE4CXgV0GcE333efhe/rpcMcdPmtBRCRX1NkBb7ezWTHwCjAY+DSE0Dl53IDVqfs7U70D/vhjnxJ26KH+AVtRUcP+AtI01AGnqaQj3nY15lbeEed19E+dK3r67IgN+3iHvGpQcgWOg7yTPWu/mQCM7vwmAHsX+n4fl/lsiAfXfAWASfMO9eO+68ft+oF3uu0/SmZNLFsJQGUyNp2LHXFzdMD1viqymbUHHgOuDCGsq/5c8BSv9V/AzM43sxlmNqOkpGTb4337wrHH+vzdp59uWPEiIpmsXh2wmRUCU4HnQgi3J4/NA44JISxNxolfDiEM3NVxao4Br18PJ5wAb73lC9uMGQOjRvnZatKy1AFniJ11xJ19dkNFT5/VsH4f72RXJWfM2YHeM31r37cBOLuzz5ronRxnXpnvN2HlEQD85YPBALR5z1dt6/qBz9Jo97Efx75YBUBIzqyr3OIdcTavvhalA06GF8YDc1Phm3gKGJP8eQzw5O6+eIcO8OyzcO21vjbDWWf5jIfzz/dhidTqYyIi2ag+syCOBF4F3gNSkXg9MB14FNgTWIRPQ1u1q2PVNgsipaICXn4ZJkyAxx7zBXS6dPFpaEcf7XOAhw6FwsLd+etJfakDzlA1O+Lk4oFVY8RJR9zfx35XD/T9yg/wNShO3ncOAKd1fgeAPfL98fe2+lUDpizzseH35/on5J3m+Of2Xeb7GX1tPvVZEpQkHXFy7bpt6xFn0RU6Ys2C+Aewsxce2VSF5Of7NLSRI+Gee+DJJ+Gll/xMuD//2fdp29ZPN06dkPGVr2i6mohkrt2aBdFYu+qAd2XZMj9DLrUOxD//6UNMZjBwoK/7MHRo1RoQuqDl7lMHnCXqGiP+knfEG/dOOuJ9fb+NA72jPaD/EgAO7+prTbTN97HduRv9Uh2vLfa1Jkrn+fzhzvP8ZTt95Ku7FS3x+cNhtXfGlZvgI99YAAAMkElEQVQ2Zc316qJ0wOmgZ0844wzfANasgddf9w/vZs3ycJ40qWr/vfaqCuShQ/2qFb176/RlEUkvGdEB18eKFR7Gqe3tt7e/OGaXLlVLUh54oN8OHqxFfFLUAWepnXXEyZl1lUlHvLmffyOs6e8fsqzv7x/3tNnL5/0Wd/Ux3gLzxz9b5x31qiU+9b/9R97LdVqYzB9e5GPB+ctWE5K5w9uu4Fyeul5dZv0Xy9kOuD66dfMpbSecUPXY+vU+XPHPf/osi3ff9fUm1q+v2qe4ePuF2YcM8cV/9GGfZIWdLAyfOpHCknBsu9Q/TGmbLOS+tZcPMazf028/2tMDd0tvP05Bx9LtblMLx5e39UjZ2skDvWOnVrRKjp23MrlMUmpReF0mKXsCuDYdOvgsiiOPrHosBFi0aMfLEz3zDCT/RykqgkGDdgxmXeFYRJpSVgdwbcy86y0u9tXWUrZuhXnzvEtOhfLf/w4TJ1bt07mzD1vUHMro2LGl/xYiDbSTjpjkRIq85NfDVsv9Yp+tP0k+bOvhHfCm3n6W1Mae/vyW1EUP2vlxy9v67cY+3qlUFBXRvr131W3b+fBHfrLIfFibDE0kV8HNlg/rdkfOBfDOtGrlgXrggds/vno1zJ69fbc8cWLVWsXgV0lOXZ4+dduzp7plEdk1BXAdunTxecdHHVX1WAjw2WfeLb/zjm9vvw1TplTt86UvVQXywQf7okMDBiiUJc2kxl2Dd52VW5LuMzmRwpITKwpW+fSyTp8lC7nv4Z3x1u6pTtg/NNnSOel8k+UEKouMTd2SBYLyk666ld8vbOUrcOWv9dttFwxNXjtTP6zbHQrgBjCDPff07bTTqh5fu9ZDedYsD+VZs+D22/3ySuDzk4cP923ECD+RpEuXOH8HEYlPAdyEOnXasVsuLYU5c2DGDJg+HaZN8/UvUj/U99vPw3j4cDj8cB8Cyav3GnUizaTmJZOSGQu20Rd+t9U+o6H1Mp/h0LpjMq2ts98v7eLjvWXt86ko9K64MrkaeWnnZO1Z85kSBYUeQ3nJ1KOQjEOHLf6NsO205ixc6EcB3MyKiqqGIb7/fX9s3ToP5GnTPJSffdanx4FPpzv2WDj+eD8tu39/DVuIZCsFcAQdO8Jxx/kGVVPjXnkFXnzRt8mT/bni4qo1Mk4+WUMWEklq9kRZMnc3GZ+11Hzitf6pdN4KXwyoTbtkvLddGyo7+IBwRZukwy1IOorkN73K9t4RW/Du2fKS5wt8dkRe6gSOLBwbVgCngepT40aP9v9X8+ZVhfFjj8H48X5yyIknwpln+rrJnTrFrlxEGkMBnIbM/ESQQYPgkkt8qc633vIgfvRR+MtffGjj5JM9jE8/XYvYSwvb2XziGrMnrKiQvDb+nzM/WSoztE7GgIuSjjh/+zG2kMyOsMoaHW7qw5FS33/bvOEMHhvWxz0ZID/fP6j71a/gk0/gjTc8mGfOhO98xxcf+tnPYOXK2JWKyO5QAGcYMw/j22+HTz+FF17wOcY/+YlfVfrSS308WaRFheBbZQVUVhDKSgllpVRu2kTlmrW+laygsmQFYVkJYVkJfLECvlhBXska8krWYOs2+ra1FNtauu3QVljoW1FqK/KtsMC3/HxfbMgs4z6xVgBnsLw8/3Du6aerLul0332w//5w550+dCEi6UsBnCUGD4bf/96X4DzmGLjqKp9X/P77sSuTnBYCobycUF5O5datvm3c5Nva9b6tXuNb0imHdRt827TJt7IyQupsJsDy83wrKMAKCnyMLumCM60TVgBnmb32gqlT4eGHYeFCD+HXXotdlYjURrMgspAZnH22L8N5/PE+de2pp3y4QiSaGutOhNTZduVJt5qaQZEsHk/e9ovJ76yrteTxYHnJ/qmXqTEGl4azJNQBZ7F+/fzkjv79/XJOS5bErkhEqlMAZ7kePeDxx32947Fj07IJkFy3kxkUodS3yi1bqdyylbA12UrLtt9CoNZLq1meb9vup9/YsAI4BwwYADff7Ff9eOml2NWISIoCOEdceKEvh3nPPbErEamnmp1xMpsilJf5VlHhW1m5b6n7lYFQ/Sy6mp1wGknPqqTJtW4NY8bAk0/Cxo2xqxERUADnlBNP9AuPTpsWuxKRRqjRGW/bQuWut5Q0GgtWAOeQww7z2zffjFuHiDjNA84hnTr5rIiPPopdiUgzyMApPuqAc0xxsS/iIyLxKYBzTLduWrZSJF0ogHNMly6wenXsKkQEdiOAzSzfzGaZ2dTk/t5mNt3MFpjZI2ZW1HxlSlNp2xY2b45dhYjA7nXAVwBzq92/BbgjhDAAWA2c15SFSfNo3VoBLJIu6hXAZtYXOBW4P7lvwHHAlGSXCcDXm6NAaVqFhVBtaVURiai+HfCdwDVAajbzHsCaEEJyVTwWA31q+0IzO9/MZpjZjJKSkkYVK41XUOAnY4hIfHUGsJmdBiwPIcxsyAuEEO4LIQwLIQzr3r17Qw4hTSgvDyor695PRJpffU7EOAL4NzP7KtAa6AjcBXQ2s4KkC+4LaLXZDKAAFkkfdXbAIYTrQgh9QwjFwFnA30II5wAvAWcku40Bnmy2KqXJmGXkCUMiWakx84CvBX5gZgvwMeHxTVOSNCcFsEj62K21IEIILwMvJ39eCHyl6UuS5pQmi0CJCDoTTkQkGgVwjlEHLJI+FMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEkm9AtjMOpvZFDP7wMzmmtlhZtbVzJ43sw+T2y7NXayISDapbwd8F/BsCGEQcBAwF/gR8GIIYV/gxeS+iIjUU50BbGadgKOB8QAhhNIQwhpgFDAh2W0C8PXmKlJEJBvVpwPeGygBHjCzWWZ2v5m1A3qEEJYm+ywDetT2xWZ2vpnNMLMZJSUlTVO1iEgWqE8AFwBDgXtDCIcAG6kx3BBCCECo7YtDCPeFEIaFEIZ17969sfWKiGSN+gTwYmBxCGF6cn8KHshfmFkvgOR2efOUKCKSneoM4BDCMuAzMxuYPDQSmAM8BYxJHhsDPNksFYqIZKmCeu53GTDRzIqAhcB38fB+1MzOAxYBZzZPiSIi2aleARxCeAcYVstTI5u2HBGR3KEz4UREIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgk9QpgM7vKzN43s9lmNsnMWpvZ3mY23cwWmNkjZlbU3MWKiGSTOgPYzPoAlwPDQgiDgXzgLOAW4I4QwgBgNXBecxYqIpJt6jsEUQC0MbMCoC2wFDgOmJI8PwH4etOXJyKSveoM4BDCEuA24FM8eNcCM4E1IYTyZLfFQJ/avt7MzjezGWY2o6SkpGmqFhHJAvUZgugCjAL2BnoD7YCT6/sCIYT7QgjDQgjDunfv3uBCRUSyTX2GII4HPg4hlIQQyoA/AUcAnZMhCYC+wJJmqlFEJCvVJ4A/BUaYWVszM2AkMAd4CTgj2WcM8GTzlCgikp3qMwY8Hf+w7W3gveRr7gOuBX5gZguAPYDxzViniEjWKah7Fwgh/BT4aY2HFwJfafKKRERyhM6EExGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRiUQBLCISiQJYRCQSBbCISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogDOMSNGwBVXxK5CRAAshNByL2ZWAixqsReU3bFXCKF77CJEckmLBrCIiFTREISISCQKYBGRSBTAIiKRKIBFRCJRAIuIRKIAFhGJRAEsIhKJAlhEJBIFsIhIJApgEZFIFMAiIpEogEVEIlEAi4hEogAWEYlEASwiEokCWEQkEgWwiEgkCmARkUgUwCIikSiARUQiUQCLiESiABYRieT/AzD1Ad6NKfzyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%% Sinkhorn\n",
    "\n",
    "lambd = 1e-3\n",
    "Gs = ot.sinkhorn(a, b, M, lambd, verbose=True)\n",
    "\n",
    "pl.figure(4, figsize=(5, 5))\n",
    "ot.plot.plot1D_mat(a, b, Gs, 'OT matrix Sinkhorn')\n",
    "\n",
    "pl.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
